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This report summarizes the result of research done under NASA 
NAG3-882 Nonlinear Mechanics of Composites with Periodic 
Microstructure . The program was carried out between 3/10/88 and 
3/8/91 with Dr. A.D. Freed of NASA Lewis acting as grant monitor. 

The effort involved the development of non-finite element 
methods to calculate local stresses around fibers in composite 
materials. The theory was developed and some promising numerical 
results were obtained. It is expected that when this approach is 
fully developed, it will provide an important tool for calculating 
local stresses and averaged constitutive behavior in composites. 
NASA currently has a major contractual effort (NAS3-26491) to bring 
the approach developed under this grant to application readiness. 

The following report has three sections. One, the general 
theory that appeared as a NASA TM, a second section that gives 
greater details about the theory connecting Greens functions and 
Fourier series approaches, and a final section shows numerical 
results . 
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Summary 

This work is concerned with modelling the mechanical deformation or constitutive behavior 
of composites comprised of a periodic microstructure under small displacement conditions at 
elevated temperature. A mesomechanics approach [1] is adopted which relates the microme- 
chanical behavior of the heterogeneous composite with its in-service macroscopic behavior., 

Two different methods, one based on a Fourier series approach and the other on a Green’s 
function approach, are used in modelling the micromechanical behavior of the composite 
material. Although the constitutive formulations are based on a micromechanical approach, 
it should be stressed that the resulting equations are volume averaged to produce overall 
“effective” constitutive relations which relate the bulk, volume averaged, stress increment to 
the bulk, volume averaged, strain increment. As such, they are macromodels which can be 
used directly in nonlinear finite element programs such as MARC, ANSYS and ABAQUS or 
in boundary element programs such as BEST3D. 

In developing the volume averaged or “effective” macromodels from the micromechanical 
models, both approaches ( i.e . Fourier series and Green’s function) will require the evalua- 
tion of volume integrals containing the spatially varying strain distributions throughout the 
composite material. By assuming that the strain distributions are spatially constant within 
each constituent phase — or within a given subvolume within each constituent phase of the 
composite material, the volume integrals can be obtained in closed form. This simplified 
micromodel can then be volume averaged to obtain an “effective macromodel suitable for 
use in the MARC, ANSYS and ABAQUS nonlinear finite element programs via user consti- 
tutive subroutines such as HYPELA and CMUSER. This “effective” macromodel can be used 
in a nonlinear finite element structural analysis to obtain the strain-temperature history at 
those points in the structure where thermomechanical cracking and damage are expected to 
occur, the so called “damage critical” points of the structure. The “exact micromechanical 
models can then be subjected to the overall “effective” strain-temperature history obtained 
at the “damage critical” location and used outside of the finite element program to evaluate 
the heterogeneous stress-strain history throughout each constituent phase of the composite 
material. This variation must be known in order to evaluate the damage history variation 
throughout each constituent phase of the composite material. 

* Work funded by NASA Grant NAG3-882. 
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1 Introduction 

The ultimate objective of this work is to produce a computer program to analyze the hetero- 
geneous stress and strain history variation at the “damage critical” locations of a composite 
structure operating at elevated temperature. This report describes some of the theoreti- 
cal foundations for the program. A mesomechanics [1] approach is adopted which relates 
the micromechanical behavior of the heterogeneous composite to its in-service macroscopic 
behavior. 

Some composites are actually comprised of a periodic microstructure whilst others are 
possessed of an essentially randomly distributed microstructure. Pictures of metal matrix 
composites (tungsten-fiber-reinforced superalloys) which exhibit a periodic microstructure 
are shown in Fig. 1 which is taken from the article by Petrasek et al [2], When the fibers in a 
composite material occupy a large volume fraction of the material, the induced deformation 
in one fiber interacts with and alters the induced deformation in the neighboring fibers. W hen 
the fibers are densely packed the interaction effect becomes dominant and must be accounted 
for in the constitutive formulation. 

At NASA-Lewis Chamis and his colleagues [3,4,5] employ two different approaches for 
analyzing the behavior of structural composites. One method employs a sophisticated finite 
element analysis of a periodic microstructure. A unit cell in the periodic microstructure is 
modelled with one hundred and ninety two three-dimensional elements and the eight nearest 
neighbor cells of the fibrous composite are modelled with superelements. By applying the 
strain-temperature history at the “damage critical” location in the composite structure to 
the superelement model, the stress-strain history throughout the unit cell can be computed 
and used to estimate the maximum damage in the composite structure. This method will 
necessarily require large resources in computer time and memory to analyze the viscoplastic 
behavior of the composite structure under in-service thermomechanical loading conditions. 

Another approach adopted by Chamis and his colleagues [3,4,5] which avoids large com- 
puter resources — is to employ composite micromechanics theory to derive simplified rela- 
tionships which describe the thermomechanical constitutive behavior of multilayered fibrous 
composites. 

When suitable boundary conditions are applied to the superelement model of the periodic 
microstructure, it is possible to predict the elastic properties of the equivalent homogenized 
material. A comparison [5] with the homogenized elastic properties predicted by the simplified 
micromechanical equations generally shows good agreement with the superelement model 
except for the Poisson ratios. At high volume fractions (~ 60%) the longitudinal Poisson s 
ratio for unidirectional fibers predicted by the simplified equations is about 15% too small, 
whilst the transverse Poisson’s ratio is about 30% too small. These anomalies occur because 
the interaction between the fibers is not accounted for in the simplified micromechanical 
model. This may be important when considering the highly nonlinear behavior of viscoplastic 
composites at elevated temperature. 

Dvorak [6] and Dvorak and Bahei-El-Din et al [7,8,9] have also made great progress in 
modelling the micromechanical behavior of nonlinear composite materials and are embarked 
on a combined experimental and theoretical effort. The variation of the stress-strain history 
throughout the unit cell of a periodic microstructure is obtained with a finite element analysis 
in which the interaction effects of the surrounding cells is accounted for by applying periodic 
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boundary conditions to the surface of the unit cell. 

Work on the theoretical foundations behind the homogenizati on of micromechanical con- 
stitutive models to produce bulk macroscopic models has been under way in Prance by Devries 
and Lene [10], Lene [11,12], Duvaut [13], Renard and Marmonier [14], Lene and Leguillon 
[15] and Sanchez-Palencia [16,17]. These references give a good account of the work be- 
ing conducted in France by other researchers. Much of this work is founded on the use of 
multivariable asymptotic techniques [18]. In an infinite periodic structure the stress-strain 
history in each unit periodic cell is, perforce, identical. Due to the finite size of the compos- 
ite structure the effects of surface tractions and displacements on the surface will cause the 
stress-strain history to vary from cell to cell. If the unit cell is much smaller than the size 
of the structure this variation from cell to cell will be small. If L is a typical dimension of 
the unit periodic cell and D is a typical dimension of the composite structure, then the ratio 
L/D is a small parameter of the problem. The displacement variation throughout the unit 
cell will then depend on six spatial variables, i.e., 

U l lij (x i , X2 , X3, X j , X % , ^3) 

where x* — x l L/D. The spatial variables x* take into account the slow variation of the 
displacement from cell to cell due to the finite size of the ratio L/D when u t is a periodic 
function of the variables x t . By expanding the displacement and other spatial variables of the 
problem into a series in powers of L/D and equating like powers in the perturbation expansion, 
it is possible to obtain the effect of the finite size of the structure on the deformation behavior 
in the unit cell. Due to the perturbative assumption of small L/D this method is not expected 
to be valid for thin composite sections or to be applicable at those places in the structure 
where surface effects or nonperiodic inclusions are important. 

Rather than employing finite element techniques to determine the stress-strain history 
variation throughout the unit periodic cell, Aboudi [19] has recently developed a macro- 
scopic formulation for periodic composites based on volume averaging a viscoplastic consti- 
tutive model over the unit periodic cell. This work expands the heterogeneous displacement 
throughout the constituent phases of the unit cell as linear and higher order functions of the 
coordinates. Good agreement with experimental results was achieved by volume averaging 
Bodner’s [20] viscoplastic constitutive model over the unit periodic cell, but the method is 
general and any constitutive model may be used to represent the deformation behavior of the 
constituent phases. The limitation here is that large spatial gradients in the strain history 
may not be accurately modelled by linear or quadratic interpolation functions on the unit 
periodic cell. 

Weng and his colleagues [21,22] have employed self-consistent methods to study the effect 
of inclusion size and volume fraction on the stress distribution in and around spheroidal 
inclusions embedded in an “effective” non-uniform matrix material, and the effect which this 
has on the overall “effective” macroscopic constitutive behavior of the composite. In the 
first paper they point out that the derivation of the fictitious body forces which represent 
the inelastic behavior of the heterogeneous composite material should be obtained from first 
principles rather than using their heuristic approach. In the second paper the Mori-Tanaka 
theorem [24] is used to represent the effect of the heterogeneous composite, and a similar 
procedure is followed in the present work to develop a self-consistent method for composites 
which exhibit a periodic microstructure. In addition, the present report also derives the 


2 


fictitious body forces for a periodic microstructure from first principles. In reference [23] 
Zhu and Weng have used a combined micromechanics and continuum theory approach to 
develop a creep deformation model for particle- strengthened metal matrix composites. They 
stress the fact that the creep resistance of the composite is underestimated when simplified 
metallurgical and mechanics approaches are adopted. 

A comprehensive application of micromechanics to mechanical deformation problems is 
given by Mura [24] in his book “Micromechanics of Defects in Solids” . This work was used by 
Nemat-Nasser and his colleagues who have exploited the mathematical simplicity of a periodic 
microstructure in order to develop elastic, plastic and creep constitutive models [25,26,27,28] 
for composite materials. The assumption of periodicity allows the heterogeneous stress, strain 
and displacement fields to be expanded in a Fourier series, which greatly simplifies the ensuing 
computations. This technique fully accounts for the interaction effects between neighboring 
fibers. Even when the composite is comprised of closely packed fibers distributed at random 
the method gives accurate results [25] for the “effective” elasticity tensor. When densely 
packed fibers form a large volume fraction of the composite material these interaction effects 
play a dominant role and must be included in the calculations. It appears that inclusion of the 
interaction effects can be as, or more, important than inclusion of the random nature of the 
microstructure when the fibers occupy a large volume fraction of the composite material. In 
this report we have developed the Fourier series approach in order to handle the viscoplastic 
behavior of the constituents in the unit periodic cell. 

The nonlinear constitutive behavior of composites with a periodic microstructure can also 
be treated with a Green’s function approach [29,30,31,32,33]. Here, the periodic heteroge- 
neous material property variation — due to the fibers — is treated as a fictitious body force in 
the matrix material. The Green’s function is used to evaluate the displacement due to a unit 
point force in the matrix material and the actual displacement at any point in the composite 
can then be determined by summing (integrating) the effect due to a volume distribution 
of fictitious periodic body forces. It is shown in Appendix B that this method is exactly 
equivalent to the Fourier series approach by invoking a mathematical technique known as the 
Poisson sum formula. The Green’s function approach is more general in that the method can 
also handle the nonperiodic case where there may be inclusions in one unit cell but not in 
the neighboring cells. It is also able to handle surface effects, although the surface integrals 
which represent the surface effects in the Green’s function method could be expanded in a 
Fourier series for thin composite sections. 

The approach adopted in the present work is to develop homogenization techniques which 
can provide simplified macromodels for use in a nonlinear finite element program, similar in 
spirit to the simplified models used at NASA-Lewis, but which account for the viscoplastic 
interaction effects in the periodic structure and which allow surface effects for thin struc- 
tures to be taken into account. Once the strain- temperature history at the “damage critical” 
location has been found from the finite element analysis, it can be used to “drive” the mi- 
cromechanical relations in order to obtain the stress-strain history variation throughout the 
unit cell. These micromechanical relations are the same relations which are used to obtain 
the simplified homogenized constitutive model. When the unit cell is chosen to have the form 
shown in Fig. 2, it is clear that a periodic arrangement of such a microstructure allows for 
the analysis of laminated composite structures. 
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2 Overview of Theoretical Modelling Approaches 

2.1 Outline of Approach 

In order to develop a homogeneous macroscopic constitutive model from micromechanical 
principles it is necessary to know the stress-strain history throughout the unit cell of a periodic 
composite. Some of the approaches which are presently being given currency are described 
in the introduction. 

The Fourier series and Green’s function approaches can be used to compute the viscoplas- 
tic stress-strain history throughout the unit cell of a periodic composite and can be simplified 
to produce a model suitable for use in a nonlinear finite element program. In this report 
the Fourier series and Green’s function approaches are developed and shown to be equivalent 
to each other by means of the Poisson sum formula. This equivalence holds only for an in- 
finitely extended medium. When the medium has a finite size the effect of spatially varying 
displacements and tractions on the surface of the medium must be accounted for. This is 
easily accomplished with the Green’s function method by retaining the appropriate surface 
integral contributions which are discarded in the case of an infinite medium. In the Fourier 
series approach the surface integral could be included, and, in the case of a thin composite 
section which has an infinite surface the integral can be expanded into a Fourier series. In 
fact, the methods can be combined, so that if inclusions are present in one unit cell and not 
in the neighboring cells, their effect can be taken into account in the Fourier series method 
by treating them with a Fourier integral or Green’s function approach. 

An overview of the present work is depicted in Fig. 3. Simplified versions of the mi- 
cromechanical constitutive equations can be volume averaged to produce a macroscopic ho- 
mogenized viscoplastic model. This can then be used in a nonlinear finite element program 
to analyze the structural behavior of a composite structure under in-service thermomechan- 
ical loading conditions. The finite element analysis yields the strain-temperature history at 
the “damage critical” location and this history can then be applied to the micromechanical 
equations to determine the heterogeneous stress-strain history throughout the unit cell. 

A detailed flow chart in Fig. 4 shows the anticipated structural analysis procedure. Both 
the Fourier series and Green s function approaches can be used to create a coarse subvolume 
model. This coarse model can then be homogenized and included in a user defined constitutive 
subroutine in a nonlinear finite element program. The Green’s function approach can also be 
used to derive a simple self-consistent model for use in the user subroutine. 

2.2 Homogenized Macroscopic Equations 

A periodic composite material is supposed acted upon by an imposed strain increment 
and responds in bulk with a stress increment Aer^. These values are then equated to the 
respective volume averaged quantities in order to obtain the “effective” constitutive relation 
for the composite material, i.e., 

= ^/// Acr o( r )^( r ) and A 4 = ^JJJ AeJ( r ) dV ( T ) (1) 

v v 

where V is the volume of the body. 
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In section 4 it is shown that the volume averaged or “effective” constitutive relation for 
the composite material can be written as 


A a% = DT iU tel, ~y /// {£>!?« Ac*, (r) - «!>„*, (r) [As£ (r) - Ac*, (r)] } dV (r) (2) 

° Vc 


where V c is the volume of a unit periodic cell in the composite material, Ae£(r) is the total 
strain increment at point r in the periodic cell due to the imposed uniform total strain 
increment Ae° kl at the surface of the composite, and Ac kl (r) is the strain increment at point 
r in the periodic cell representing the deviation from isothermal elastic behavior. The fourth 
rank tensor 6Dij k i(r) is defined by the relation 


SD ijkl ( r) = m (D{ jkl - D™ kl ) 


(3) 


where $(r) = 1 in the fiber and i?(r) = 0 in the matrix, with D f ijkl denoting the elasticity 

tensor of the fiber and D™ ki that of the matrix. 

In the expression for the average or “effective” constitutive relation in equation 2, the 
quantities Ae° w D™ kl and 8D ljkl (r) are given. The deviation strain increment Ac fc; (r) can 
be obtained throughout the periodic cell as a function of position r by using an explicit 
forward difference method since the stress and state variables in a viscoplastic formulation 
will be known functions of position at the beginning of the increment. Everything is therefore 
known explicitly except the total strain increment Ae^(r). 


2.3 Fourier Equation Overview 

In the Fourier series approach described in section 4 we find that the total strain increment 
is determined by solving the integral equation, 

1 ±°o 

A eh (r) = Ae° kl + tt 5^ 9kiij (C) x 

V C Tlp=0 

X /// e * Ac ^ (r ' } “ 6D *” (r ' } [ Ae « (r ' } _ ACrs (r,) l } dV (r,) (4) 

v c 


where the fourth rank tensor g k uj (C) is given by 

gmj «) = l (<j (C) + CAM*' (0) 


(5) 


in which the Christoffel stiffness tensor M {j (C), with inverse Af y 1 ( 0 is defined (cf. [33]) by 
the relation, 

M l j (0 = D™ jq Cp(q (6) 

with Cp = ip/y/^m^m = ip/i being a unit vector in the direction of the Fourier wave vector £, 
and £ = v // £m£m denoting the magnitude of the vector £. In equation 4 the sum is taken over 
integer values in which 

27rni . 27rn 2 . 2 t m 3 (7 \ 

Cl — t ) S2 t ) S3 j ' ' 

L 1 L, 2 *-'3 
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and Li, L 2 , L 3 are the dimensions of the unit periodic cell in the X\, x 2 , x 3 directions, so that 
V c = LiL 2 L 3 . The values of n u n 2 , n 3 are given by 

n p = 0, ±1, ±2, ±3, . . . , etc., for p = 1,2,3 (8) 

and the prime on the triple summation signs indicates that the term with ni = n 2 = n 3 = 0 
is excluded from the sum. 

2.4 Green’s Equation Overview 

In the Green’s function approach the total strain increment Ae^(r) is determined by solving 
a different integral equation, viz., 


A£y(r) — ^ £ fcl + JJJ Uklmn ( r r ) ( F ) 

V 

- fiDmnrs (r') [AeJ s (r') - A c rs (r')] } dV { r') 


(9) 


where the fourth rank tensor U k imn ( r — T ') gives the hi component of the total strain incre- 
ment at point r due to the ran component of a stress increment applied at point r' in the 
infinite matrix with elasticity tensor D™ nrs , i.e., 


Ukimn (r r ) 


1 ( d 2 G tm (r - r') + a 2 G, m (r-r') \ 

2 \ dxidx n dx k dx n ) 


( 10 ) 


and the volume integration in equation 9 extends over all the periodic cells in the composite 
material, i.e., over the entire composite. 

The Green’s function tensor is defined in Appendix A, equation A. 26, by the Fourier 
integral [24,32,33] 

d 3 K Mjj 1 (C) e _jK.(r-r') 


Gij < r - r<) = III 


(27t) 3 


--1 

~K* 


( 11 ) 


in which the tensor £ is now defined by the relation Q = Ki/ K with K — yj K q h q denoting 

the magnitude of the vector K = {K\,K 2 , K 3 ). 

In Appendix B it is shown, by applying the Poisson sum formula, that equations 4 and 9 
are identical, although the summation extends over the integer values n u n 2 , n 3 in equation 4 
and extends over the periodic cells in equation 9. 


2.5 Integration of the Equations 

Both equations 4 and 9 are implicit integral equations for the determination of the total strain 
increment Ae^)(r), since this unknown quantity appears both on the left hand sides of the 
equations and on the right hand sides under the volume integrations. 

The “effective” constitutive relation given in equation 2 and the total strain increment 
relation, given by either equation 4 or 9, contain the volume integration of the deviation 
strain increment Ac^;(r). In the periodic cell the deviation strain increment at any point r 
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will be determined from a unified viscoplastic constitutive relation [34] appropriate to the 
constituent phase in which the point r resides. If a constituent phase is included at the fiber- 
matrix interface, a constitutive relation can also be proposed for this phase, and the resulting 
inelastic strain increment determined for inclusion in the volume integrals. This may be 
important for metal matrix composites where there can be chemical reactions between fiber 
and matrix at elevated temperatures, and for composites where the fiber has been coated to 
enhance overall composite properties. 

Equations 2, 4 and 9 form the basic incremental constitutive equations for determining the 
“effective” overall deformation behavior of a composite material with a periodic microstruc- 
ture. In order to update the stress state in each of the constituent phases in preparation 
for integrating the “effective” constitutive relation over the next increment, the constitutive 
relation, 

A <Tij( T ) = D m { r) (Ae£(r) - Ac«(r)) (12) 

is used, where Aj«( r ) = D {jki or ^ijki according as the point r is in the fiber or matrix. 
This relation is used to update the stress <7,.,(r) and, in turn, the internal viscoplastic state 
variables q , ( r ) at each point r in preparation for computing Acfc/(r) in the next increment. 

The derivation of the preceding equations and some methods for their solution are dis- 
cussed in the succeeding sections of this report. Numerical solutions will be obtained during 
the research effort from appropriate FORTRAN computer programs. 

3 Periodic Microstructure 

3.1 Volume Averaging 

The periodic composite is supposed acted upon at its surface by a spatially linear displacement 
increment, Ati°(r), given by 

A u°i (r) = Xj Ae°. + xjA uj% (13) 

where Ac*) and Acc*) are the spatially uniform strain and rotation increments at the surface 
of the composite. 

If the matrix material was homogeneous and had no fibers embedded in it, the strain 
increment would be homogeneous and given by 



Since this is constant, we may trivially volume average Ae° over the volume V of the homo- 
geneous matrix material to obtain 

_ 1 rtn ( 9 ( Au °( t )) s ( a “?M) 

u VJJJ 2 ( dx s + dXi 

which, by Gauss’ divergence theorem, may be written as 

A 4 = ^ JJ l ( n l( r ) A 4( r ) + f*i(r)AuJ(r)) dS( r) (16) 

s 
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where the integral extends over the surface of the material and n,( r) denotes the outwardly 
directed unit normal vector at point r on the surface. Thus, by a pply ing the displacement 
increment Au°(r) in equation 13 over the surface of the material to produce the surface strain 
increment given in equation 16, equations 15 and 16 show that the strain increment in the 
matrix material is spatially uniform. 

If the displacement increment Au°(r) in equation 13 is applied to the actual composite 
material, the total displacement increment within the material, A uf (r), will vary in a periodic 
manner due to the assumed geometric periodicity of the composite material, so that 

Auf(r) = Au°(r) + Atq(r) (17) 

where Auf(r) is the displacement increment which would be induced in the homogeneous 
matrix if the fiber phase were absent, and Aufir) is the perturbation or deviation from the 
homogeneous value due to the presence of the fibers. 

Corresponding to these displacement increments, the total strain increment at any point 
r in the composite, Aeffir), is given by the relation 

Ac£(r) = AfiJ, + Ae w (r) (18) 


where 


A e° kl = 


'd{A u° k ) + d(A uf)' 


dxi 


dx k 


and 


A _ l(d(Au k ) d(A Ul y 
Aew(r) = 2 


(19) 


with AejJ/ representing the spatially constant total strain increment which would be produced 
on the surface and in the interior of the homogeneous matrix if the fibers were absent, and 
Asfci(r) representing the deviation from the uniform value due to the presence of the fibers. 
Both the total strain increment Asffir) and the perturbed strain increment Ae w (r) vary 
throughout the composite in a periodic manner. 

We define the volume averaged stress and strain increments as Acr^ and Ae® , respectively. 
The required “effective” constitutive equation for the composite material is then an expression 
relating the volume averaged stress and strain increments. For a function /( r) which varies 
with position the volume average is defined by the relation 

(f) = ^JJJf(r)dV(r ) (20) 

V 


Since the composite is 
we may write 


assumed to be comprised of a periodic aggregate of identical unit cells, 


(/) 


1 

V c 


JJJf(r)dV(r) 

v c 


( 21 ) 


where V c denotes the volume of the unit periodic cell. 

If we volume average the total strain increment in equation 18, we obtain 


(A el) = f Uf Ael(r)JV(r) = AeJ, + f J [ f A £ „(r) dV (r) (22) 

c v; c v c 
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<Ae£) = AeJ, + (A**,) (23) 

But the volume averaged total strain increment is defined as Ae*,, so that (Ae*,) = Af*, and 

(As*,) = 0 (24) 

which shows that the volume averaged perturbation strain increment, Ae*,(r), is equal to 
zero. 

3.2 Eigenstrain and Deviation Strain Increments 

If the elasticity tensor is denoted by Dij k i( r) and the inelastic strain tensor by e kl ( r), then 
the constitutive equation at any point r in the composite material can be written as 

0ij( r) = D ijkl { r) (4i(r) - e£(r) - a k i{r) ( T - T 0 )) (25) 

where a k i( r) is the coefficient of thermal expansion. 

The incremental form of Hooke’s law is 

Aa tj (r) = Dijkit r) (Ae£(r) - Ac«(r)) (26) 

where Ac^r) denotes the incremental strain representing the deviation from isothermal elas- 
tic conditions and is given by 

Ac kl (r) = Ae^r) + a* kl (r)AT (27) 


in which 


a* kl (r)AT = a k i{r)AT + (T-T 0 )Aa k i(r)- 

- D k ilj(r) AD ijmn (r) (s^Jr) - e£ n (r) - a mn { r) (T - T 0 j) (28) 

is the nonisothermal increment in strain. The tensors AD ijk i(r ) and Aa kt (r) represent the 
incremental changes in the elasticity and thermal expansion tensors due to the temperature 
increment AT. 

In a unified viscoplastic constitutive formulation [34] which is integrated by an explicit 
forward difference method, the inelastic strain increment Aejy(r) is a function of the current 
stress (at the beginning of the increment), ^(r), and the current values of the internal 
viscoplastic state variables, c/ t (r). For example, if 

£ij = fij Q 3 ) ( 29 ) 

then Ae[j = fa (cr rs ,q 3 ) At, and the inelastic strain increment is independent of the total 
strain increment Ae^j(r). This independence of the inelastic strain increment on the total 
strain increment is no longer true if an implicit integration method ( e.g . backward difference) 
or subincrementation method is used. 

The elasticity tensor D ijkl (r) may be written as 

£>„«(■ r) = CJ*, + «Aj U (r) (30) 
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where 


«W r) = m (D' jk , - D? Jkl ) (31) 

and $(r) = 1 in the fiber and $(r) = 0 in the matrix, the superscripts / and m referring to 
the elasticity tensor of the fiber and matrix, respectively. The constitutive equation at any 
point r can then be written, from equation 26, as 

A<7ij( r ) = (D™ kl + 6D ijU (rj) (Aeg, + Ae kl (r) - Ac fcZ (r)) (32) 

or 

Aa y (r) = D? m (Ae° + Ae«(r)) - 
- {£>"„ Ac„(r) - «Wr) [AeJ, + Ae„(r) - Ac„(r)] } (33) 

If the quantity in braces is set equal to D^ kl As k[ (r), that is, if 

D iJkA £ ki( r ) = D ijkA c ki{r) ~ 6D ijkl ( r) [Aeg, + Ae w (r) - Ac w (r)] (34) 

then equation 26 can be written in the form, 

Ao-y(r) - D™ kl (Ae£(r) - Aefc(r)) = D™ kl (Ae°, + Ae*,(r) - Ae^(r)) (35) 


From the preceding equation it is evident that the eigenstrain increment, Ae^(r), represents 
the incremental deviation from isothermal elastic behavior in the composite material when 
the elasticity, tensor is taken to be a spatially constant tensor appropriate to that of the 
matrix phase. 

Newton’s law for continuing static equilibrium throughout the strain increment requires 
that 

dAa^r)) 


dxj 


0 


Equations 35 and 36 then require that 


(36) 


or, if Ae k[ is constant, 


dxj 


™ d(Ae M (T)) _ 

U ijkl q x U ijkl 


Ae;,(r))} Q 

(37) 

dAeUr)) 

dxj 

(38) 


4 Fourier Series Approach 

4.1 Fourier Expansions 

The application of Fourier series to the calculation of the “effective” overall constitutive 
behavior of periodic composites has been dealt with in detail by Nemat-Nasser and his col- 
leagues [25,26,27,28]. This work is used in this section to develop constitutive relationships 
for viscoplast.ic composite materials under small displacement conditions. 
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Due to the geometric periodicity of the composite we may expand Au k (r) and Ae kl {r) in 
a Fourier series (c/., for example, Appendix 3 of Mura’s book, [24]). This gives 

+ „ j-oo +oo /2irni 27 rn 3 \ 

Au fc (r)=£ EE' Ia " l 1 *) (39) 

m=0 ri2=0 ri3=0 

where L u L 2 , L 3 are the dimensions of a unit cell in the x u x 2 , x 3 directions. The coefficients 
Au k in the Fourier expansion are determined by multiplying each side of equation 39 by 

/2777711 27T7712 2777713 A 

e - l {~LT Xl+ ~LT X2+ L 3 X3 J anc j integrating over the volume of the unit cell to give 


Ll L 2 £3 / 277711 

Au fc (n 1 ,n 2 ,n 3 ) = / / j Au k {v)e 

t » =0 t o D rro=f) 


<277711 277712 277713 


+ -lT X2+ ~lT X3 


) dx 1 


dx 2 dx 3 (40) 


x 1 =0 £2=0 £3=0 


where only the terms with m, = n» survive in the summations. 
Equations 39 and 40 can be written in shortened form as 

±00 

Au t (r) = £££ AMO**' 

n p =0 

with coefficients A u k (£) determined by the inverse relation 


A u k (£) = ^ IIJ Au k {r)e * r dV( r) 


where 

^ = (6,6,6). T = {xux 2 ,x 3 ), V c = LiL 2 L 3 (4dj 

with 

£• = - (no sum on i ) for i = 1, 2, 3. (44) 

The strain increment Ae kl (r) can also be expanded in a Fourier series to give 

±oo 

(45) 

7l p =0 

with coefficients A £& determined by the inverse relation 

Ai kl (£) = yJJJ A4(r)e- i4r dV (r) (46) 

c v c 

In equations 41 and 45 the prime indicates that the term with rq = n 2 = n 3 = 0 is excluded 
from the summations, since Au k (ni = 0, n 2 = 0, n 3 = 0) represents a rigid body displacement 
increment and Ae* kl (m = 0,n 2 = 0 ,n 3 = 0) represents a spatially uniform strain increment. 
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4.2 Equilibrium Equation 

By substituting equation 41 into equation 19; equation 19 into the left hand side of equation 
38; and equation 45 into the right hand side of equation 38, the equilibrium relationship 
becomes 

±00 1 / v 

D T,m E E E' o ( Afl * « > + A “< «) ) e i{ r = 

n p = 0 Z V 7 

±oo 

-^EEE' A ^(5)^ iir 

Tlp—O 

or 

(4) = -iD? jk &te* kl (0 (47) 

If 4 = VMi denotes the magnitude of the vector £, a unit vector £ in the direction of £ can 
be written as = 4i/4- Equation 47 can therefore be written in the form, 

eDT ia (|) (!) A * t «) = «) 

or 

e (A"«00) A “fc «) = -iBTwtAK , «) (48) 

The second rank tensor, 

M* (0 = Uk, (C) = D.’JhOC, (49) 

is called the Christoffel stiffness tensor (c/. [33]) and equation 48 can be written as 

?M lh (C) «) = -«£>£.£, Ai (€) (50) 

This equation can be inverted by premultiplying each side by the inverse tensor £ -2 M -1 to 
give the Fourier expansion coefficients 

An, «) = -M,l' (<) D? lr £AK, (0 r 2 (51) 

The expansion coefficients can now be substituted into the Fourier expansion of Au^r) in 
equation 41 to give 

±oo 

Au k (r) = - E E E' ' (0 Air. «) r (52) 

n p =0 

This result may now be substituted into equation 19, so that the perturbation strain increment 
may be written as 

±00 1 / \ 

Ae«(r) = E E E O (C) (A + C 2 M,T' (<) fcft ) D” .Air*. (€) e is r (53) 

n p =0 Z V ' 

If we define the fourth rank tensor g k uj (4) by the relation 

9m (0 = \ (m* (<) (jO + M,r' (0 6&) (54) 
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( 55 ) 


then the perturbation strain increment can be written in the form 


Ae„(r) = E E E' m (0 D^M:. «) r 

n p — 0 

and by inserting the relation for the Fourier expansion coefficients A i* 3 from equation 46, we 
obtain 

A EH (r) = ^EEEWC) [f[D™,te‘ r Ar')e , « r -' , ' l dV(T') (56) 

Vc n p =0 Ve 

where the integration extends over the volume, V c = L1L2L3, of the unit periodic cell. 

From equation 18 the total strain increment is given by 



(0 jJjD^AKA^e'^-^dVir') 

V C 


(57) 


which, from the definition of D% r , AeJ, (r') in equation 34, may be written in the final form, 


A4M = A4 + iEEE'tod(C)//A ,t<r - r 'HA>.Ac rI (r')- 

- 6D tjr , (r') [a 4 (r-j - Ac., (r')] } dV( r') (58) 

This implicit integral equation — equation 4 in section 2.1 — must be solved to yield the total 
strain increment Ae^(r) at each point r in the unit periodic cell. 

Instead of solving for Aejy(r) from this implicit integral equation, we could use equation 34 
to eliminate Ae^(r) from equation 57 to give an equivalent integral equation for Ae^(r), viz . , 

0” w A4,(r) = -D”„Ac t ,(r) - «D««(r) [AeJ, - Ac„(r)] - 

-HWr)±E£E's*""><<) ///D™„AC(r')c‘«<- r ')<fV'(r') (59) 

c n p =0 y c 

The incremental constitutive relation at any point r is given in equation 35, and this 
relation can be used to update the stress state at any point r in the unit cell once equation 59 
is solved for Ae^(r). Alternatively, equation 58 can be solved for Aer^(r) and inserted into 
equations 34 and 35. The overall “effective” constitutive relation for the composite material 
can be obtained by averaging equation 35 over the unit periodic cell. This gives 

(A<r ij ) = (l>5 u (A4 + Ae H -A4)) 


or 

(Aery) = £>”« A4 + D? iU (Ae fc i) - (A 4) (60) 

If we equate the volume averaged stress increment (Aery,) and the overall bulk response stress 
increment A<j° , he., if (Aoy,) = A<r t °, and we note from equation 24 that the volume averaged 
perturbation strain increment is zero, i.e. (A £ki) = 0, then the overall “effective” constitutive 
relationship is 

A ff » = DT m AcJ, - D? jkl (Aei) (61 ) 


13 


or 




i -yjjj {D? lt A c u (r) - 6D m (r) [A el, (r) - Ac*, (r)] } dV( r) (62) 


which is the result presented in equation 2 of section 2.2. The procedure for integrating the 
overall “effective” constitutive relation then proceeds as follows. 


4.3 Fourier Integration Algorithm 


1. From a knowledge of the stress state throughout the unit periodic cell at 
the current time, t, calculate the inelastic strain increment Ae£ (a rs , q s ,r) from an 
appropriate unified viscoplastic constitutive relation. The viscoplastic constitutive 
relation will vary according as r is in the fiber or matrix phase, respectively. 

2. Compute the eigenstrain Ae* kl (r) throughout the unit periodic cell from the 
implicit integral equation 59 or from equations 34 and 58. 

3. Compute the stress increment throughout the unit periodic cell from equa- 
tion 35 and update the stress, strain and viscoplastic state variables according to 
the relations 

a v} (r, t + At) = a l} (r, t) + Aa l3 (r), 

4 (r, t + At) = ej (r, t ) + AeJ(r), 
q t (r, t + At) = q x (r, t) + Aq t (r) . 

4. Calculate the overall “effective” stress and strain increment for the compos- 
ite from equation 61 and update the overall “effective” stress and strain from the 
relations 

4 (t + At) - 4 (t) + a < t °, 

4 ^ = 4 (0 + A 4 - 

5. Repeat the preceding calculations for each incremental load step. 


4.4 Implicit Integration Algorithm 

The preceding algorithm makes use of the fact that the inelastic strain increment Ae kl (r) is 
independent of the total strain increment Af4(r) if an explicit forward difference method 
such as Euler or Heun forward difference — is used to integrate the unified viscoplastic relations 
for the fiber and matrix phases. If an implicit method — such as backward difference or sub- 
incrementation — is used, the inelastic strain increment depends on the total strain increment. 
In this case the total strain increment must be obtained by iterating equation 58 in the form, 


A4)(r) = A4, + ^Y.T.T:9mAOjlJ^ r -’^D^,Ac,,( r \Ael,(^ 

c n p = 0 v c 

- SDij rs (r') [AeJ s (r') - Ac rs (r', Ae^ (r'))] } dV(r') 


(63) 
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The first iterative guess can be taken as Ae^(r) = Ae^, and the right hand side evaluated 
to give an improved guess for Ae^(r). This process is then continued with 


{ A£ «M} j+1 


+ 77 EEE' am (0 y//‘«* £(r - r ' > {DZrAcr. (r', { AeJ, (r')}J - 

C n p ~ 0 y c 

- 6D ijr , (r') [{a£ (r')} A - Ac„ (r', {AeJ, (r')}J] } dV(r') (64) 


until the A th and (A + l) th iterates of Ae^r) converge. 

Equation 59 is not so convenient for iteration as equation 58 when Aef s (r) depends on 
Ae^(r). It is always necessary to know the total strain increment Ae^(r) in order to calculate 
the inelastic strain increment Aef s (r', Ae^ (r')). But equation 34, viz., 


D™ kl A4i(r) = D? ]kl Ac kl (r, Ae T pq {v)) - 6D tjkl (r) [Ae£(r) - Ac fci (r, Ae£(r))] (65) 


is an implicit equation for Ae£,(r) when the iterated quantity, Ae* kl (r), is given. Equation 63 
is therefore the appropriate equation to iterate when the inelastic strain increment depends 
on the total strain increment. 

The procedure for solving the implicit integral equations in 58, 59 and 63 is described in 
section 8. 


5 Green’s Function Approach 

5.1 Green’s Solution of Navier’s Equilibrium Equation 

The equation of continuing static equilibrium for the composite material throughout an ap- 
plied strain increment is given by 

d(A* tJ ( r ) ) + = Q (66) 

dxj 


where A/,(r) is the incremental body force per unit volume of the composite material. From 
equations 35 and 66 we obtain 


n m 

U ijkl 


d(Ae&( r)) _ d 


dxj 


dXj K Ae « (r) ) 


A/»(r) 


(67) 


which is equivalent to equation 37 in the absence of the incremental body force Afi(r). Prom 
this equation it is clear that the divergence of the stress variation produced by Ae^(r) may be 
formally regarded as a fictitious body force increment, analogous to A/*(r), which is applied 
to the homogeneous matrix material with elasticity tensor D™ ki . The theory of elasticity for 
homogeneous materials is generally concerned with the solution of the homogeneous differ- 
ential equation 67 — Navier’s equation — when the right hand side is zero. When body forces 
are present the standard method of solution is to obtain the displacement solution at r due 
to a unit body force applied at r'. This solution is given by the Green’s function Gij (r - r') 
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which gives the displacement in the z th direction at r due to a unit point force applied in 
the j th direction at r'. For a distributed incremental body force A/* (r') the displacement 
increment at r is obtained by summing the results for the distribution in the form 


A«i(r) = JJJg y (r - r')A /, (r') dV( r') (68) 

V 

The integration extends over the whole volume, V , of the composite material which may be 
regarded as being of infinite extent. 

When A/j (r') = 0 we know that the displacement solution is Auf(r) = Au°(r), cor- 
responding to an applied uniform strain increment Ag:F on the infinite boundary of the 
homogeneous matrix. For an effective distributed body force increment given by the right 
hand side of equation 67, with A f i (r') = 0, the solution for the total displacement increment 
A uf (r) can be written as 

A uj (r) = Au°(r) - JJjG ik (r - r ')^y (^l„A4 n (r')) dV ( r ' ) ( 69 ) 

v 1 

This corresponds to equations 17 and 39, the volume integral corresponding to the perturbed 
displacement increment Aiq(r) in 17. 

For a material which is homogeneous with elasticity tensor D™ kl the Green’s function 
satisfies the differential relation (cf. Appendix A, equation A. 11), 

^ + ^0-r')=° P») 

w r here 8 im is the Kronecker delta tensor given by 8 im = 1 if i = m and 6 im — 0 if i ± m, and 
8 ( r - r') is the three dimensional Dirac delta function defined by the relation 

6 (r - r') = 6 (aq - x[) 6 (x 2 - x' 2 ) 6 (x 3 - x' 3 ) (71) 


By applying the Fourier integral techniques in Appendix A, the Green’s tensor is shown to 
have the Fourier integral form, 


Gij (r 



d 3 K Mij 1 (C) 

(2t r) 3 K 2 


g-jK.(r-r') 


(72) 


in which the inverse Christoffel stiffness tensor (cf. [33]) M tJ 1 (C) is defined by 

M ,- 1 (C) = feCpC ,)" 1 ( 73 > 

with C P = K P lyjK m K m = K p jK being a unit vector in the direction of the Fourier wave 

vector K, and K = J K m K m denoting the magnitude of the wave vector K. 

Making use of the relation 

Gik (r — (j^Zmn^ £ * mn (r (r — r )D klmn Ae mn (r )^ — 

- aGik f x 7 r,> PS„„AC„ (r') (74) 
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we may write equation 69 in the form 

AuJ (r) = A«J(r) - /// A (<?* (r - (r')) dV(r') + 

V 1 

+ dV( r ') (75) 

V 1 

The first volume integral can be transformed into a surface integral via Gauss’ divergence 
theorem, viz., 

Ill lb ( G ‘‘- (r " r ') D S- A£ ”" < r '0 = II n ‘ ( r ') G -‘ ( r " r ')DS ra „Ae^ (r') dS( r') 

(76) 

The surface integral extends over the entire outer surface of the “infinite” matrix material. 
Since this is assumed to be at an infinite distance, all the integration points r' in the surface 
integral are at an infinite distance from the field point r and Gik (r — r') = 0. Thus, for an 
infinite body the first volume integral in equation 75 vanishes. This would not be the case 
for a finite body in which the field point r is close to the surface integration point r' , and 
the volume (or surface) integral would need to be retained for these situations. In this case 
other surface integrals would arise (c/. Appendix D, equation D.27) due to the application of 
boundary incremental displacements or surface tractions on the surface of the material. 
From the properties of the Green’s function, 

dGjk (r - r') _ dG ifc (r - r') 
dx[ dxi 

which follows since G & is a function of 

r - r' = (xi - x\ , x 2 - x ' 2 , x 3 - x' 3 ) 

Equation 75 may then be written alternatively as 

Auf(r) = Au?(r) - jfJ ^^^- D^ mn Ae' mn (r') dV( r>) (79) 

V 1 

But AeJ(r) = \{d (Auf(r)) /dxj + d (AuJ(r)) /dx t ) , so that by differentiating equation 79 
with respect to x x and Xj and taking half the sum, we obtain 


(77) 

(78) 


AeJ(r) = A4 + JJJu ijU (r - r') Ae' m „ (r') dV( r') 

V' 

which, by means of equation 34, may be written as 

AeJ(r) = A e% + JJJ U m (r - r') {D£„ Ac,, (r') - 

V 

- 6D Ur , (r') [a£ (r') - Ac,, (r')] } dV( r’) 


(80) 


(81) 
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An equivalent integral equation, involving the eigenstrain increment Ae* (r), can also be 
obtained by using equation 34 to eliminate Ae^(r) from equation 80, which gives 


DTjkA^ii r) = Ac w (r) - 6D ilkl ( r) [AeJ, - Ac w (r)] - 

- SD.M r) JJJu, (r - r') D^AC, M dV(r') 

V 


In the preceding equations the operator, 

1 / a 2 G*(r-r') | &G jk (r-r f ) 

2 y dxjdxi dx.idxi 



(82) 


(83) 


gives the ij component of the strain increment at point r due to an applied stress increment 
component kl at point r' in an infinite homogeneous medium with elasticity tensor D"j k , and 
Green’s function given by equation 72. 


5.2 Equivalence of Perturbed Strain Increment 

From equations 18, 56 and 80 we see that the perturbed strain increment, Ae w (r) = Ae^(r) - 
Ae'h , is given by the equivalent relations, 

1 ±OC 

V C n p ^ 0 

or 

Ae k ,(r) = JfJ U Um „ (r - t>) £>”„» Ac), (r') dV(r') (85) 

The volume integral in the Fourier series representation extends over the volume, V c , of 
the unit periodic cell and the summation extends over the integers n p = 0, ±1, ±2, . . . , etc., 
where p — 1,2,3. In the Green’s function approach the volume integral extends over the 
entire infinite medium, he., over all the periodic cells comprising the material. It is shown in 
Appendices B and C that the Fourier summation expression in equation 84 can be converted 
into the Green’s function expression in equation 85 by means of the Poisson sum formula. 

From equation 34 it is evident that if the elastic properties of the fiber are the same as 
that of the matrix, then 6D tJ ki( r) = i9(r) (D{j k i — Dijki) — 0> * n which case 

A4(r) = Ac fci (r) (86) 

is known explicitly without having to solve the integral equation. From equations 58 and 81 it 
can also be observed that Ae^(r) is known explicitly when 8D m { r) = 0. The explicit relation 
in equation 86 holds only when an explicit forward difference method is used to integrate the 
viscoplastic constitutive relations. For implicit integration methods in w^hich the inelastic 
strain increment Ae^(r) depends on the total strain increment Ae^(r), equations 58 and 
81 show that even when 6D iik i{ r) = 0, the equation to determine Ae£(r) is still an implicit, 
integral equation. 


jJJ D2 nr ,AC, (r')c‘ s "- r ' 1 dV{ r') (84) 
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6 Self-Consistent Method 


6.1 Outline of Self-Consistent Method 

In this section we establish a self-consistent relationship between the overall “effective’ stress 
increment, Act", and the applied strain increment, Ae° , for a matrix material which ha s 
cylindrical fibers embedded in it in a periodic fashion. 

From equations 34 and 61, this relationship can be written as 


AaJ = D™ kl Ae° kl -y JJJ {D% kl Ac k i (r) - SD ijkl (r) [Aef, (r) - A c kl (r)] } dV( r) (87) 

c Vc 

where the total strain increment is determined from equation 81 in the form, 


Ae£(r) = Ae° kl + jjj U klrnn (r - r') {^ nrs Ac rs (r') - 

V 

- SD m „, M [a£ (r-) - Ac„ (r')J } dV( r') 


( 88 ) 


These equations can be solved in an approximate fashion by means of a self-consistent 
method in the following manner. 

First, assume that the unit periodic cell consisting of a cylindrical fiber embedded in a unit 
matrix cell, Fig. 5, is replaced by a cylindrical fiber (of radius = a) embedded in a cylindrical 
matrix (of “effective” radius = b) as depicted in Fig. 6. The other unit cells outside the given 
unit cell — -i.e., the rest of the composite — are then smeared out into a uniform matrix material 
whose overall “effective” constitutive properties are the volume average of the constitutive 
properties of the constrained unit periodic cell. The “effective” constitutive properties will 
be transversely isotropic if the fibers are arranged in hexagonal arrays or tetragonal if they 
are arranged in square arrays. 

Second, assume that the total strain increment, Ae^(r), and the strain increment repre- 
senting the deviation from isothermal elastic behavior, Ac k i(r), are spatially constant in the 
fiber and matrix phases of the unit cell. These constant values (different in the fiber and 
matrix of the unit cell) are taken to be the volume averages over the respective constituent 
volumes of the fiber and matrix phases of the unit cell. 

The composite now consists of three constituent phases, viz., the fiber, matrix, and 
smeared out average phases. If the elasticity tensors of these phases are denoted by D,j k i i 
D™ kl and D ijkl , respectively, then the elasticity tensor at any point r in the composite can be 


written as 

D k lrs(r) = Dklrs T bD k i rs { r) 

(89) 

where 

bDklrsi. r) ^^klrs ^ klrs ^klrs 

(90) 

if r is in the fiber; 

bD k lrs( r ) = fiD™l T s = Dklrs ~ ^ klrs 

(91) 

if r is in the matrix; and 

6D k i rs (r) = 0 

(92) 
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if r is in the surrounding smeared out “effective” material. 

The fiber and matrix constituent phases now represent fictitious body forces in the infinite 
“effective” medium with elasticity tensor T>ij k i, and the total strain increment is obtained from 
the solution of the integral equation, 

AeJ.(r) = Ae% + JJJu ijkl (r - r')D klrs Ae* rs (r') dV( r 7 ) (93) 

v 

in which 

D klrs Ae* rs (r 7 ) = D klrs Ac rs (r') - 6D klrs (r 7 ) [a£ (r 7 ) - Ac rs (r 7 )] (94) 

6.2 Strain Increments in Three Phases 

We now make the approximation that the strain increments in the three phases are spatially 
constant and equal to their respective volume averages, so that if r' is in the fiber A ejj (r') 
and A Cij (r') are replaced by 

A ej(/) = ~ JJJaeJ, (r') dV( r') (95) 

f Vf 

and 

A c, (D = y jj! Ac,, (r') SV(v') (90) 

1 Vf 

so that, from equation 27, 

Ac„(/) = AeZ(f) + a-'AT (97) 

with 

a,! AT = a{,AT + (T - To) Aa{, - 

- AD' lmn (eUD - £,(/) - <*L V - To)) (98) 

If r 7 is in the matrix the relations arc replaced by 

AeT (m) = ~ JIJ Ael (r') dV(r') (99) 

m 

and 

Ac,j (m) = Tjjj Ac., (r') dV( r') (100) 

m v m 

whore 

Acij(ra) = A ef/m) + a*™ AT (101) 

and 

a*™ AT = a™ AT + (T — T 0 ) Aa| — 

- (A™*/) 1 AD klmn (4mM - £ mn( m ) ~ “1 ( T “ T o)) ( 102 ) 
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If r' is in the smeared out “effective” or homogenized medium the corresponding results are 



A^ - Ae% = TF HI A^l ( r ') dV ( r ') 

3 V s 

(103) 

and 

15 lj = yJJjAc„(r')dV(r') 

s V, 

(104) 

where 

A Cij = A + a-jAT 

(105) 

and 

a* AT = oiijAT + (T - T 0 ) Aa {j - 



— Dijkl A D k l mn — £mn ~ a rnn (T ~ Tq) j 

(106) 


The volumes Vf, V m and V s refer to the volumes of the fiber, matrix and smeared out medium, 
respectively. If V c 1® the volume of the unit cell and V denotes the total volume of the entire 
composite material, then 

V c = Vf + V m and V = V c + V S = Vf + V m + V s (107) 

6.3 Applied, Homogenized and Volume Averaged Increments 

At this point it is important to emphasize the following distinctions. First, the strain incre- 
ment applied to the composite is denoted by Ae° which causes an incremental stress response 
To obtain the overall “effective” constitutive equation these are equated to the cor- 
responding volume averaged quantities, ^A and (A c^). In the “effective homogenized 
medium all quantities are denoted with overbars. 

At any point r the appropriate constitutive relation is 

Acrjj(r) = D ijk i{T) (Ae£(r) - Ac w (r)) (108) 

If we volume average this relation over the unit cell we obtain 

(A Gif) = (D ijkl Aeh) - (D tjkl Ac kl ) (109) 

In the homogenized phase the constitutive relation can be written as 

Ac jj = Dij k iAe kl Dij k i Ac k i (HO) 

Since the strain increment Ae k( in the homogenized phase must correspond to t he a pplied 
strain increment As^— as in equation 103— and the homogenized stress increment Ac tJ must 
correspond to the overall bulk stress increment A c?, we write the constitutive relation for 

the homogenized phase as 

Ac° = D ijkl A4 - D tjkl Ac kl (111) 
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6.4 Requirement for Self-Consistency 

For self-consistency we require that the volume average of the microscopic constitutive relation 
in equation 108 over the unit cell, viz. equation 109, should correspond to the constitutive 
relation for the overall “effective” homogenized medium in equation 111. That is, 

^°ij ~ ~ (DijkiAcu) = DijkiAeh — D ljkl Ac k i (112) 

for self-consistency. Under the approximation that the strain increments A £ kl and A c kl are 
spatially constant in the constituent phases, we obtain 


Acr° = y JJJ Dijki(r) (Ae£(r) - Ac«(r)) dV( r) = D ija i 


As kl — DukiAc , 


ijkl L 


~kl 


or 


(113) 


= y&UH ( A 4(/) ~ *cM)) + ( A 4i(™) - Ac t ,(m)) 

= DijkiAe° kl — DijkiAc k i (114) 

At this point the elasticity tensor D 1}k i and the deviation strain increment A c k i in the homoge- 
nized medium are unknown quantities. In the next section we will solve the integral equations 
for the total strain increments in the fiber and matrix phases, A s kl (f) and Ae kl (m), and we 
will find that these values depend on the quantities Dqfcj, Ae kl and A c k i in the surrounding 
homogenized medium. Then, by equating the coefficients of Ae kl on both sides of equa- 
tion 114 we obtain a relationship for the unknown elasticity tensor D ijk i of the “effective” 
homogenized medium. The value of the unknown deviation strain increment A c k i in the ho- 
mogenized medium can then be obtained by equating the terms independent of Aejj on the 
left hand side of equation 114 to the corresponding term Dij k iAc k i on the right hand side. 

We now obtain the total strain increments Ae kl (f) and A £ k i(m) in the two phases of the 
unit cell. First, consider the total strain increment in the fiber phase. 


6.5 Total Strain Increment in Fiber Phase 


Equation 93 can be volume averaged over the fiber phase to give 

yJJJ 5(r)dV(r ) = A4 + i JJJdV( r) jffu m (r - r') D u „Ae% (r') dV( r') (115) 

f V S f Vf V 

where the field points r are in the fiber volume, Vf, and the integration points r' are in all 
three volume phases (V = Vf + Vm + V,). Equation 115 can be written as 


H(f) = A 4 + v, 


/// dV(r) /// UiiU (r “ r<) (/) dV >') + 
. V, { Vf 

+ III u ijkl (r - r') D klrs Ae* rs (m) dV( r') + JJJ U m (r - r') D klTS Ac rs dV( r') 


(116) 
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in which D klrs Ae* s (r') has been replaced by 

Dkirs Ae* s (/) = D k i rs Ac rs (f ) - SD[ lrs 

A e T „(f) - Ac„(/)] 

(117) 

and 



DkirsAers (m) = DkirsAcrsim) - 6D™ lrs 

A ej s (m) - A c rs (m) 

(118) 

in the respective fiber and matrix phases, and by 



& kirs^Crs 

(119) 


in the smeared out “effective” medium where SD k i r3 (i - ') = 0. 

In the first integral in equation 116 the field point r lies in the volume Vf, and since 


JIJ U ijkl (r - r') dV(r') D klrs = S ijrs (120) 

v 

is Eshelby’s tensor (cf. Appendix E, equation E.l and [35]), which is a constant tensor 
independent of r when the field point r lies within the cylindrical volume V included in an 
infinite medium with elasticity tensor D ktrs , we may write the first integral as 

UJu ijkl (r - r') dV( r') D klrs Ae* rs (f) = S ljrs Ae* rs (f ) (121) 

Vf 

The second volume integral extends over the volume V m = V c — Vf of the matrix phase. Thus, 
for the second integral, 

III Uljkl ^ ~ dV ( T> ) Dkirs A£* rs (m ) 

V m 

= JffUijki (r - r') dV(r') D klrs Ae r * s (m) - JJJu ijk i (r - r') dV( r') D klrs Ae* rs (m) 

v r Vf 

= Sj jrs AE* s (rn) - S ijrs Ae* r ,(m) = 0 (122) 

since the field point r lies in the cylindrical volume Vf and therefore within the cylindrical 
volume V c . 

We now have to deal with the last integral in equation 116. This integral can be written 

as 

JJjUijkt (r - r') dV{ r ) D k [ rs Ac rs 

v s 

- _ fff 1 ( d dGuJjj-V) d dGjkjjj- V) 

.//./ 2 \dxj dxi &Xi dxi 

Vs 

[[[l ( - ^fc( r ~ r/ ) , d dGjkjr-r') 

JJJ 2 \dxj dxi dx^ dxi 

V .s 


dV (r) D k i rs Ac rs 


j dV( r') D k i rs Ac rs 
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which can be transformed via Gauss’ divergence theorem into two surface integrals: one over 
the outer “infinite” surface S of the smeared out “effective” medium; and the other over the 
inner surface of the “effective” medium, i.e., over the surface S c of the unit periodic cell. The 
volume integral then takes the form, 


JJJ U ijU (r - r 7 ) dV (r') D ktrs Ac rs 

v s 

- Ill 
+S fll 

Sc 

where n* (r 7 ) is a unit normal vector at point r 7 on the surfaces pointing away from the volume 
V s . Now since the field point r lies in the fiber and the integration points r 7 on the surface S 
are infinitely removed, we have dG ik (r - r')/dxi — » 0 on the outer surface S of the composite, 
and the first surface integral can be neglected. If we write n* (r 7 ) = — n* (r 7 ), then n, (r 7 ) is a 
unit normal vector pointing away from the volume V c on the surface S c of the unit cell, and 
we have, via Gauss’ divergence theorem and equations 77 and 83, 


V) ”“£ - ) -+"f(O ac<t(r ~ ,,) 


n *j c n 


dxi 

, N dG lk (r - r') 


dxi 


+ < (r') 


dxi i 

t\ 9Gjk (r — r 7 ) ' 


dxi 


dS{r ) D klrs Ac rs + 
dS(r') D k [ rs Ac rs 


IJJ u -*‘ ' r - r') dV(r') D klrs Ac, 

= ~H\ ("^ <r '> Bx, 


A 0G lk (r r 7 ] + n _ (r , } OG jk (r r' ) ^ 


dxi 


- /// 

V c 


1 ( d dG, t (r - r') , d d G Jt (r - r') 


2 \ d xd 


dx. 


+ 


dx'i dxi 


dV(r') D kirs Ac r 


1 /a 2 G lt (r-r') , d 2 Gjk (r - r') 


+ 


III 2 V dxjdxi ' dxidxi 

= - JJJ Uiju (r - r 7 ) dV (r 7 ) D ktrs Ac TS 


dV{r') D klrs Ac r 


Vc 


Since the field point r lies in the cylindrical volume V c , the preceding equation takes the form, 


or 


JJJ Uiju (r - r 7 ) dV (r 7 ) D klrs Ac rs = ~JJJ U ljld (r - r 7 ) dV (r 7 ) D klrs Ac r 

V', Vc 

JJJ U ijk i (r - r 7 ) dV{ r 7 ) D k i rs Acr S = -S ljrs £c rs 


(123) 


where S ljr!i is Eshelby’s tensor for a cylinder with elasticity tensor D ljk i. 
From equations 116, 121, 122 and 123 we obtain 

AsJjif) = Ae®j + — JJJ {Sij rs Ae rs (f) — <9jj rs Ac r .,| dV(r) 

f v f 


I 


1 
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and since the quantity within the braces is a constant tensor, 

A ej(/) = Ae% + S ijr , (A - S^.) 

Now 

DurAK, U) - Dki r ,ACr,(f) ~ ? lD llr, [AsJ.U) - Ac„(/)] 


(124) 


vso that 

Ae;. (/) = A c„(i) - D,^ 6D [AeK/) - Ac ra „(/)] 

and on substitution into equation 124 we obtain, 

A 4 (f) = Ae“ + S,,„ (a Cr.U) - Ai r , - D^SD^ [A 4 „(/) - A c_()) 

Given that 

^ijkl — 2 “I" 

denotes the fourth rank identity tensor for symmetric second rank tensors, the* preceding 
equation can be written as 

\jijmn + ^ijrsD rspg SD^ H ?mn ] A £ mn (f) = ^ e ij + 

+ S{j rs (a Crs(f) - Ary s ) + SjjrsD rspq ^DpqTnn^ C Tnn(f) 
which, by premultiplying each side with the inverse of the tensor in square brackets, gives 


~ lijmn d~ Sij rs D rspq 6Dp qmn J |Ae mn + Sm nr3 ^Ac r5 (/) Ac rs ^ + 
~k SmnklDklpq SD^ACr.U)} 


(125) 


The phase volume averaged stress increment in the fiber is then given by the relation 


A<r , ,-(/) = D’ ijkl (A - A c H (/)) 


(126) 


6.6 Total Strain Increment in Matrix Phase 

Now consider a field point r in the matrix phase. From equation 93 we may write 

AfJ(r) = Ar° + jjj U, jU (r - r') dV(r') D tlr ,Ae‘ r ,(f) + 

Vj 

+ JIf U ' jkl (r “ r ') dV(r')D klrs Ae* rs (m)+ Jjju ijkl (r - r') dV{v')D klr ^ rs (127) 

v. 

Since = V r — V/, the second integral can be written as 

JIJ U ijkl (r - r') dV{ t') D klrs Ae; s (m) 

v m 

= jjj u ijkl (r - r') dV( r') D klrs Ae* ra (m ) - JJJ U ljkl (r - r') dV{r') D klrs Ae* rs (m ) (128) 

i;. v, 
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and from equation 123 the last integral in equation 127 may be written as 

JJJ Uijki (r - r') dV (r') D klrs Acr S = -JJJ U ijk i (r - r') dV (r') D kWs Ac rs ( 129 ) 

v, v c 

so that equation 127 is transformed into 

A £ J(r) = Ae° + JJJu ijkl (r - r') dV(c')D Ur , (A £ ;.(/) - A e‘„(m)) + 

Vf 

+ JJJ u ijkl (r - r') dV(v') D klrs (A e* r 8 (m) - A^ s ) ( 130 ) 

Vc 

Averaging equation 130 over the matrix phase gives 

A sTjM = AeJ + y-jfjm r) ( JIJUw,, (r - r') dV(v')D klr . (Ae‘ r ,(f) - A<(m)) 


+ 




+ JJJ u ijkl (r - r') dK(r') D Wrs (A<£(m) - A c rs ) 

V, 

or, since V m = V c — V}, 

A ej,(m) = A £ » + i 


(131) 


U 1 y 

y m 


JJJ dV (r) JJJ Uij k i (r - r') dV (r) D Hra (A e* rs {f) - Ae* s (m)) 

v c v f 

+ III dV ^ JJJ Ui i kl ( r ” r ’) dV(r)D klrs (a ,e* a (m) - A c^) - 

v v c 

- fjfdV(T) HI U m (r-r')dV(rj D u „ (&€*„(/) - A e' r ,(m)) - 

'Vf Vf 

-JJJ d.V( r) JJJ U ijk i (r - r') dV(r') D ktrs (A £* s (m) - Ac™) 


+ 


(132) 


Now consider the first integral in equation 132. We may interchange the order of the 
volume integrals so that 

JJJdV(r)JJJ U iik i (r - r') dV( r') = JJJ dV( r') JJ f U m (r - r') dV (r) (133) 

i; v } vy ^ 

Now r and r' are dummy integration variables, so that on the right hand side of equation 133 
the variables may be replaced with the integration variables x and y, viz., 

JJJ dV (r) JJJ Ui jk i (r - r') dV(r') = JJJ dV( y) JJJu ijk i (x - y) dV(x) (134) 
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But, from equation F.5 of Appendix F, 

Utjki (x - y) = Uij k i (y - x) 

so that 

jjj dV (r) jjj U ijk i (r - r') dV{ r') = jjj dV(y) JJJu ijk i (y - x) dV{x) (135) 

V c Vf Vf Vc 

and the dummy variables y and x may be replaced with the variables r and r' to give 
jjjdVMjjju.Mr r') dV{ r') = r') dV(r') (136) 


v, 


Vf 


This relationship is discussed in Mura’s book ([24], page 336) where it appears under the 
heading of the Tanaka-Mori theorem. 

From equation 136, the first integral in equation 132 is integrated over the field points r 
within the cylindrical volume Vf. Since these field points lie within the cylindrical integration 
volume V c . the first integral in equation 132 may be written as 


V, 


-jjj dV( r) jjj U ijkl (r - r') dV(r') D klrs ( A e* ra (f) - A e* 5 (m)) 
” Vf v c 

y jjj S ijrs (A e*„U) - Ae* s (m)) dV( r) 


_Yi 

V m 


Stjr, (A e’,(f) - A 


( 137 ) 


In the second integral in equation 132 the field points r lie within the cylindrical volume V c 
and so the second integral may be written as 

±jjj dV(r) jjj U l]kl (r - r') dV( r') D klrs (A e* rs (m) - Ac^) 

m v c v c 

= y- jjj Sijrs (Ae* a (m) - Ac rs ) dV( r) 


= Sij rs (A e* rs (m) - A c rs ) 
v w 

In the third integral the field points r lie within the cylindrical volume Vj and so 

y jjj dV( r) jjjUijki (r - r') dV(r’)D klrs (A e* r ,{f) - A^(m)) 

' m 'Vf v> 

= y jjj Sijrs (A e* rs (f) - Ae* s (m)) dV( r) 


( 138 ) 


= T7~ Sijrs (A £*„{/) - A 

v m 


( 139 ) 
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Finally, in the fourth volume integral, the field points r lie in the cylindrical volume V}. 
Since this lies within the cylindrical integration volume V c , we have 

y JJJ dV{ r) HI U ijkl (r - r') dV(r’) D klrs (A e r * s (m) - A^,) 

m V; V c 

= y III Sijrs (Ae* s (m) ~ A c™) dV( r) 

= (A £ ;,(m) - SE.,) (140) 

* m 

We thus obtain from equations 132, and 137 to 140, 

A efjim) = Ae^ + yS ijrs (A e* a (f) - A e* a (m)) + 

* m 

+ Sij rs ^A£ r5 (?7l) Ac^^ 

* m 

- (Ae* s (/) - Ae* a (m)) - 

- §ijrs (A < a (m) - AC,,) 

^ m 

= Ae° + Vc ~ — S^s (Ae* s (m) - Ac™) 

* m 

or 

A ejj{m) = Ae^ + Sy™ (Ae* s (m) - Ac™) (141) 

This relation for the total strain increment in the matrix phase is similar to that for the fiber 
phase given in equation 124. By following the steps leading to equation 125 the expression 
for Ae^(m) can be put in the form 

Aefj(m) — j/jjrmj + Sij rs D rspq 6D™ mn J |Ae^ m + S mnrs (A c™(ra) — Ac™) + 

4“ S mnk lD fcipq «D^,Acv,(m)} (142) 

The phase volume average stress increment in the matrix is then given by the relation 

A<r y (m) = Dg u (A e^(m) - Ac fc/ (m)) (143) 

6.7 Overall “Effective” Constitutive Relation 

As stated in section 6.4, for self-consistency we require that the volume average of the con- 
strained micromechanical constitutive relation over the unit periodic cell should correspond 
to that for the “effective” homogenized medium. From equation 114 we require that 

AffJ = h A / w (a el(f) - A .<*(/)) + ^ (A el(m) - A c„(m)) 

= Dij k iAe kl - D ijk iAc k i (144) 
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where the total strain increments in the fiber and matrix phases are given by equations 125 
and 142 as 


1 M-4 



~ ^ Sij rs D rsp(] 5D^ qrnl ^ |Ae mn + S mnrs (^Ac rtt (f) A7’r.v) + 




+ SmnklD klpq ^I^pqrs^ C ra (/) j 

(145) 

and 

A efj(m) = 

lijmn d" Sij rs D rS pq + S mnrs ^Ac r5 (?7l) — Ac n j + 

+ S mn kiD klpq 6D pqrs Ac rs (m ) j 

(146) 

with the deviation strain increments defined in equations 97, 101, 105 as 




Acy(/) = AeJ(/) + a*/ AT 

(147) 



Acjj(m) = AeJ(jji) + a‘”Ar 

(148) 

and 


Acij = Ai^ + a* AT 

(149) 


By inserting equations 90, 91, 145, 146 147, 148 and 149 into equation 144 the relationship 
for self-consistency requires that 


Acr°. = Ai jk iA£ 0 kl - {AijkiSkirs Ae^ s - B ijra Ae^ s {f) - C ijrs A£^ s (m)} - 
— ^AijkiSklrs^rs ~ Bi 3TS Q* T { — Cij rs a*™^ AT 

(150) 


(151) 


(152) 


(153) 


in which 


DijiciAsfci D t j k iAe kl T , 3 k \ ol AT 


A, jki t,- Djj rs j^/rsw + S rsp< jD pqmn (p mnkl -D fnn ^.j^J + 


V, 


v,„ 


+ ~^ D i] rs Irski + S rspq D pqmn (p"' lnkl ~ Amifcf)] 


and 


B lt ki — 


y ^ijrs Irski 4 ” S rspqB pqrnn (^D mn!f/l J X 


X 


Sghkl I SghnlJI a brd (j-^cdkl Talk: l ) ] A Dl 


v c ijkl 


Cijkl = 


V, 


771 T\m 
Y u ijrs 


mngh &mngh ^ j 


-1 


Irski + SrspqDpy^ ^ D , 

Sghkl + S g j ia hD abc(i ^ D c dkl ~ Hcdkl )] y-Dijkl 


X 


V m 
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(154) 


These results for B l j k i and Cijki can be simplified somewhat. We write B^ki as 


_ V f 

Bijkl = -yVijkl 


where 

and 

We then obtain 


Vijkl — D fjpq ( Ipqrs 4 Xpq TS ) ( S rs kl 4 X rak l) 


ijkl 


Y — c n 1 sn f 

s^pqrs — ^pqmn^ mngh UJ ^ ghrs 


Vijkl + ijkl ~ &ijpq {Ipqra + Xpqrs) (Srskl + Xrskl) 

° r -1 

(Ajh) (vklmn + ^Ln) = (A/« + X ijkl) 1 (S k lmn 4 4mn) 
This result simplifies to 

( D Ijki) m mn + I ijmn — (-ft ijkl 4" Xijkl) {Sklmn 4 X k lmn) 
which can be premultiplied by {Ipqij + Xpqij) to give 

{Ipqij 4 Xpqij) ijkl) Vklmn Sp 


Jpqmn Ipqmn 


from which 


Vijkl — Dijpq ( Ipqrs + Xpq TS ) ( S rsk l Irskl) 

From this result we find that B ljk i and C ijk i can be written in the simplified forms 


Bijkl — y B ijpq ^ f pqgh 4 Spq rs D rsmn ^ C T1 ~ in gf i C rnTli jii ) J {Sghkl Ighkl) 


and 


Cijki — —D^pq \lpqgh + Spq rs D rsTnn (^D™ ngh Omnj/i)] {Sghkl Ighkl ) 


(155) 

(156) 

(157) 

(158) 

(159) 

(160) 
(161) 

(162) 

(163) 


Equating the coefficients of Ae° kl in equation 150 for self-consistency then requires that 

Dijki = Aijki (164) 

which, from equation 151, produces the implicit relation 


D{jkl — -,t Bjjrs [frsM 4 S rsp qD pq mn {j-^mnkl Hmnkl ^ j 4 

+ Yn l£> m 

' y ijrs 


Iraki *"i“ S rS pqDpq mn m nkl 


(165) 


The value of homogenized “effective” elasticity tensor Dijki may be obtained from this im- 
plicit relationship by iteration. Naturally, when the self-consistent method is embedded in a 
nonlinear finite element program, this iteration would be done outside of the code and the 
explicit values of D ljk i would be used in the program. 
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(106) 


For self-consistency we also require from equation 150 that 

Aijpq Spqfd Af kl ~ ~ k i(m) — Dijkl&Efci 


and _ 

AijpqSpqki&£i — BijkiAaJ — Cijki&OL kl — Dij k ia k[ ( 167 ) 

which, by setting A t j pq — D tJpq , reduce to 


and 


— &ijpq {^pqrs J-pqrs )] (B rM &£u(f) + Crskl^£kl( m )) 

&ij = \j-^ i ]pq ( Spq r <t Ipqrs ) j (^^rskl^kl ^ ^rskl^kl ) 


(168) 

(169) 


The overall “effective” constitutive relation for the homogenized composite in equation 1 50 
is now easily computed. 

If a forward difference algorithm is used to evaluate the viscoplastic strain increments, 
the only implicit equation which occurs in the formulation is that for the elasticity tensor 
of the homogenized medium given in equation 165. It is, perhaps, ironic that in deriving 
the highly nonlinear viscoplastic constitutive relationship for the homogenized medium, the 
only iterative procedure required is that for the elasticity tensor. This implicit elasticity 
relationship also occurs in the subvolume method due to the occurrence of the tensor 6D t jki 
in the volume integration. The implicit nature of Dijki is due to the fact that the homogenized 
elasticity tensor is found by volume averaging the constrained elastic properties of the unit 
periodic cell, and these constrained properties, in turn, depend on the elasticity tensor D,jki 
of the homogenized constraining medium. 

The constitutive relations given in equations 126 and 143 are used to update the stress- 
strain history in the constituent phases, whilst equation 150 is used to update the stress-strain 
history in the homogenized self-consistent medium. These relations, which contain AeJj(f), 
A z^j(iv) and Dijki, depend on the Eshelby tensor Sij rs for the homogeneous smeared out 
medium, which is defined in equation 123 as 


Sljrs = III Uijkl (r “ r ' } dV{ - T,) Dklrs ( 1 70) 

v 

when the field point r lies within the cylindrical volume, V . The “effective” homogeneous 
smeared out medium for a composite with cylindrical fibers will exhibit transverse isotropy if 
the fibers are arranged in hexagonal arrays, and it is shown in Appendix E that the Eshelby 
tensor for a transversely isotropic cylinder, whose cylindrical axis is normal to the plane 
of transverse isotropy, has the component form, 
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(171) 

(172) 

(173) 



S 2233 

_ -D1133 

2 Z?im 

( 174 ) 

£ll 33 

= $2233 

( 175 ) 

*^2211 

— $1122 

( 176 ) 

*S'l212 

a 3 T >1111 — Du22 

- *1221 - RTi 

oL> mi 

( 177 ) 

S 2323 

= <$2332 = <Sl 313 = S133I = 4 

( 178 ) 


with all other S t jki = 0 . If the fibers are arranged in tetragonal arrays, the Eshelby tensor 
will exhibit tetragonal symmetry. This case is currently being worked out. 


7 Integration of Self-Consistent Model 

Fourth rank tensors can be written in Voigt notation as matrices and second rank tensors as 
vectors, (cf. Appendix 2 of Mura’s book, [ 24 ]). For example, with the notational changes, 

A(7i = Ao - )!, A 02 = A022) A03 = A033, A04 = Ao'23, Ac^ = Ao^, A<T6 = A012 


and 


A^i = Ash, A £“2 = AS 22 j Ae 3 = Aff 33 , Ae4 = 2Ae23) ASs — 2 AS 13 , ^6 — 2 A£i 2 


Hooke’s law for an isotropic elastic medium can be written as 
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For a transversely isotropic medium — such as the smeared out “effective” matrix for hexagonal 
fiber arrays— the relationship can be written as 
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from which the isotropic results can be recovered by_ taking D nn = 2 /ii/(l -v)/(l - 2 v), 
D3333 = 2pt/(l - v)/{l - 2 u), Dm2 = - 2 u), D n33 = 2 \ivj{\ - 2i/), and £>1313 = //. 
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The Eshelby tensor [35] relates the constrained strain increment, in an inclusion 

which undergoes a transformation or eigenstrain increment, A e k i, in an infinite medium with 
elasticity tensor 2 in the form 

Ae-j- = SijktAeh ( 181 ) 

In Voigt notation we have 

A< = S lJ Ae* 

where the Eshelby matrix takes the form, 
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*(1 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 


G) 


0 


0 

0 

0 

0 

3£>mi — D i\22 


SDu 


u 


The integration of the self-consistent model then proceeds as follows. 


(182) 


1. Initialize the starting variables: time t = 0; temperature T = T 0 ; overall 
‘‘effective” stress and strain erf = ef = ef = 0 for i = 1 to 6; stress and strain in the 
respective phases cr l (f) = ej(f) = ef(/) = 0, and (Ji(m ) = ej (m) = ef (m) = 0 
for i = 1 to 6; equilibrium stress in the respective phases 0,(/) = f \{m) = 0 for 
j = 1 to 6; drag stress in the respective phases K(f ) = K 0 (f) and K{m) = K 0 (rn). 

2. Compute the overall “effective” elasticity matrix iteratively from the relation 

D.J = K + (dL/ - T>-i )] + 

+ Y e D * K + SuV ‘’'‘ ( DZ > - ' 

where S kj is the Kronecker delta matrix, 6 kj = 1 for k = j and 6 k j = 0 for k ± j, 
and the Eshelby matrix S tJ is given in equation 182. 

3. Start the loading history step. Evaluate the inelastic strain and state vari- 
able increments in the fiber and matrix phases from the unified viscoplastic con- 
stitutive relations. Any unified viscoplastic model may be used. Such relations 
may have the following form: 
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In the fiber the inelastic strain increment is 

A ef (/) = 


( s,(f) - n,(f)\ 

yjl (§«.(/)-«.(/)) (§*,(/) - n.(/)) 

\ K{f) ) 

K(f) 


71^ — 1 


where the equilibrium or back stress increment is calculated from the relation 

A «,(/) = efAef(/) - e ' v 'fA E f(/)A<(/)fi i (/) 

for j = 1 to 6, and the drag stress increment from 

A K(S) = - ei ( K(f ) - K„(f))} v'lAeft/JAeff/) 

In the matrix a similar set of constitutive relations can be used, so that 


Aef (m) = 


{ Si(m ) - 

y| (KM - <W)) (KM - fi <?M) 

V *("0 ; 

K(m) 


n m -l 


AQi(m) = tfAef (m) - g^y/fAef (m)Aef (m) fti(m) 
for ?' = 1 to 6, and 

AK{m) = - gj 1 (K(m) - tfo(m))] v /§As| > (m)Aef (m) 

The quantities n^, n m , gf,, g™, for p — 1 to 4, are material constants associated 
with the unified viscoplastic constitutive relations. 

The deviatoric stress in the fiber phase is defined by 

Si(f) = 0i(f) — I + a 2{f) + vzifi) f° r i= 1,2,3 

and 

Si(/) = for i = 4, 5, 6. 

Similar relations apply to the matrix phase, viz., 

Si(m) = cr l {m) — | (cri(m) + cr 2 (m) + £T 3 (m)) for i — 1,2,3 

and 

Si(m) = ai(m) for f = 4, 5,6. 

4. Compute the “effective” inelastic and thermal strain increments from the 
relations . _ 

Ai/ 5 = [D ip (Sn - /„)] “ (B qk Ae^(f) + C qk Ae p k (m)) 


and 

where 

and 


a* AT = [D w (Sn - /„)] 1 (B qk a* k f + C qk a* k m ) AT 


B q k — ~rjrD qi [<5 ip + S tm D mn (Dh p T) np jJ (*5 p fc <5 p fc) 


K + SimD-l ( D ™ - D np )] _1 (S pfc - V) 


\ 1 — 1 


V c 
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pm 


r>. Evaluate the deviation strain increments from the relations 


Ac,(/) = A <(/) + a‘/ AT 
Acg(m) = Ae^(m) + a* m AT 

n«<l P 

A c q = Ae q + a* AT 

6. Evaluate the phase volume averaged total strain increments 

AeJ(f) = [6i, + S lp D-'(D' j -T>„)}-'{Ae° + 

+ S„ (A c,U) - Ai,) + S jp D^' {p{r - Dr) A <*(/)} 

and 

As’(m) = [i>, J + S, p O„ 1 (n™-5, i )]" 1 {A e ” + 

+ Sj, (A c,,(m) - So,) + s„,73 m ‘ (d” - 75,, ) Ac r (m)} 

7. Calculate the stress increments in the fiber and matrix phases from the 
relations 

Aa,(f) = D! J (AsJ(f)-Ac,(f)) 

and 

Aa t {m) = D™ (A ej(ra) - Acj(ra)) 

8. Compute the overall “effective” stress increment from the relation 

Ao” = D„ (Ae» - Ai- - a" AT) 

9. Update the variables: 

a,(f,t + At) = <r i (f,t) + A<7 i (f) : 

ai(m,t + At) = a l (m,t) + Aai(m) 

q, (/, t + At) = a(/,*) + Aa(/) 

Qi(m,t + At) = £li (m, t) + AQ t (m) 

K(f,t + At) = K(f,t) + AK(f ) 

A' (m, t + At) = A' (m, t) + AK(m) 
e?(fJ. + At) = sf(/,t) + A £ f(/) 
sf (/77, t + At) = ef (m, <) + Aef (rn) 
e[{f,t + At) = ej (/, t) + Aef (/) 
ff (m, t + At) = ef (m, £) + Aef (m) 

^(t + At) = of(0 + Ao? 
e°(* + At) = e?(«) + Ae? 

T(t + At) = T(t) + AT 

10. Start new load step. 
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8 Subvolume Method 

8.1 Approximate Integration of Integral Equations 

The determination of the stress and strain increments throughout the composite material 
requires the solution of the integral equations 


or 


A"ti A4( r ) = D m Acm(r) - 6D ijk ,( r) [Ae^, - Ac u (r)] - 

- r)i Y. £ £'»<".» (0 (r')e‘«'- T ’>dV(r') 

VC n P=° J Vc 

^7jki^ £ ki( T ) = D™ k iA c id(r) — f>Dijki(r ) [Aej^ — Acw(r)j — 

- 6D ijU ( r) fjju tlm „ (r - r') DZ„AK, M W<J) 


(183) 


(184) 


at each field point r in the unit periodic cell. 

Nemat-Nasser and his colleagues [25,26,27,28] have demonstrated the efficacy of dividing 
the unit cell into a number of subvolumes and assuming that A4„ (r') is replaced by 

(185) 


AC„(r') = Ae± = ~ ///A £ ;„,(r') dv ( r ') 

0 V 0 


corresponding to its average value in the /3 th subvolume. 

Let there be N subvolumes in the unit cell, with M subvolumes in the fiber and N — M 
subvolumes in the matrix. Then the preceding integral equations can be written as 

Dijki^ £ ki( T ) = ^ijki^ c ki( r ) ~ fiDijkifo) [Ae fcl — Acfe/(r)j — 

DZ„. Aef. (186) 

V c rip-0 ' - I 


?=' \ Vp 


and 


D? jU AeJ,(f) = D% k ,Ac t ,(r ) - 6D ijU ( r) [AeJ, - Ac„(r)] - 
- 6D ijU (T)YY (fjj Ju*™ (r - r') rfV(r')J DZ„ r . ,Ae% 


(187) 


where V^ q denotes the / 3 th subvolume in the 9 th unit periodic cell, and it is assumed that the 
field point r is in the first periodic cell for which q = 1 . 

These equations can be volume averaged over the a th subvolume in the unit periodic cell 

to give 

D J r, t A4? = Km Act, - SD?m [Ae” u - Ac£] - 6D% U X 

{vJJJ e<(rdV{r) ) % v> {vJII e ' itr,dVir,) ) (188) 
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and 


D^,Aetr = OJ? t( A<S, - 6 Df jkl [M - Ac?,] 

iD m E E [vJII dV(l) III Ukt "» < r “ r ') dV '< r ') 


.Ad? 


£> m 

mnrs r5 


(189) 


In these equations the deviation strain increments Ac£, are evaluated from the unified 
viscoplastic constitutive relation for the a th subvolume based on the stress value ofAf) or 
a'j'jlm) in the subvolume, according as the a th subvolume is in the fiber or the matrix phase, 
respectively. The notation 6D? jkl also denotes the value of D{ jkl - D™ kl or 0 according as the 
a th subvolume is in the fiber or matrix phase, respectively. 

If we use Nemat-Nasser’s notation and write 


Q a (i) = ~Jjje^dV(r) 

Vet 


and denote 


fQ= V« 
J Vc 


(190) 


(191) 


as the volume fraction of the a fch subvolume, then the preceding equations may be written as 


N 

E 

0=1 


±oo 


DT lr .6° e + sDfa e E E'»«m« (o Dz^fcr (o q*(-o 


n v =0 


= D$ u Ac£ - 6D? jU [Ac°, - Ac?,] 


Ad? 


(192) 


and 


v 

E 

A 1 


^^tyii^iii Ukimn (r r ) dV (r ) j D rnJI 

q ~ V V <» V 0q / 

DTj,A4: - 6D U [Ml - Ac?,] 


A e*J! 


(193) 


where i = 1 if a = /? and 8 nii = 0 if a ^ (3, and no sums on a, ft are intended unless 
explicitly stated. 

Now 8D°j kl = 0 if the a th subvolume resides in the matrix. In this case equations 188 and 
189 show that 

AeJ, ft = Ac£ for M < a < N (194) 

Thus, only M unknowns (associated with the subvolumes in the fiber) are involved in Ae*f 
and the N — M known quantities (associated with the subvolumes in the matrix) given by 
equation 194 may be taken over to the right hand side of the equations. Equations 192 and 
193 may therefore be written in the form 
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M 

E 

0 = i 


±oo 


A?,.*”' 5 + «)g H E E E W (O DZ^fQ* (OQ* (-«) 

Tip —0 

= A”«Ac£ - <5A% [AeJ, - Ac?,] - 

JV ±oo 

- E *A5«EEEW(OA^./'<r(«Q' , (-«)A<& 


Ae* s 


(195) 


and 


M 

E 

jj=l 


/ mnrs 


Ae*J 


DTir.i ’- 6 + «A“« E \yI!! dv{T) III Uu ™ (r “ r,) dnt,) D; 

9=1 \ “ V a V 0v j 

= DT jU Ac?, - 6D° jk , [A4 - Ac?,] - 

- E 6D°uf;\YfffdV(r)[f[u Um „(r-r')dV(T'))DZ„ r ,A4, (196) 

0^M+l <?=1 \ Q vL / 


for a = 1 to M. 

By defining the fourth rank tensor A ( ^ ra and the second rank tensor 6? in the Fourier 
series representation as 

±oo 

A?fr. = D? jn 6 °f> + lDf m E E E'sh™ (0 DZ, r .fQ° (0 Q” (-«) ( 197 ) 

Tip — 0 

and 

6” = D^Aci, - «5g« [Ae2, - Ac?,] - 

N dhoo 

E ■5A%EEE'»<™(0^„„/‘’0“«)l? , (-eAc?, (198) 


N 

E 

/3=M-hl 

or, in the Green’s function representation, as 


and 


<, = + 6 D° u e(,' jVki™ (r - r') dV(r') (199) 

9=1 \ Q V 0 V 0Q J 


b% = DZjuAcl, - 6D° jkl [Ac?, - Ac?)] - 

- E 6D? j k,Y\Irfff<n r (r)jjJ u t , m „(r-r')dV(r')\DZ nr .,Ac?. (200) 

fl=M+l 9=1 V Q V„ 4, 
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th(' integral (equations can be cast in the form, 


M 


Y A rfrs^rs = b ?j for <* = 1 to M 

0=1 


(201) 


In Voigt notation the fourth rank tensor A^f rs can be written as a matrix A'^f, and the 
second rank tensors Ae*f and bf, can be written as vectors Ae*^ and b " , so that 


M 


Y A^Ae*P — b^ for a = 1 to M 

0=1 


( 202 ) 


This represents a system of 6M linear equations for the unknown values A e*^, each matrix 
element a/3 of the matrix A consisting of a 6 x 6 submatrix, in the form 


[An] 

[A 21 ] 

[Aja/] 

• • [a 2M ] 


{Ae* 1 } ' 
{Ae* 2 } 


{b'l 1 
{b 2 } 

' • ■ ■ [A 0 / 3 ] 


< 

{Ae*^} 

> = < 

{b"| ’ 

[Aa/i] 

■ • • [Amac] 


, K"} , 


. M . 


(203) 


where the submatrix elements are defined as 

Ag Ag Ag A g A 


-A-a/?] 


A a P 

A n 

a<*P 

-* 1 21 

A Q & 

^31 

A 


afi -| 
16 


aol(3 A a @ A a P A a P 4 a/? 

^22 ^23 ^24 ^25 71 26 


ac*{ 3 A a & A a P A a P A a & 

Aw ^33 /I34 ^35 ^36 


l 32 

A a & 

"M2 


A a(3 A 

/I43 


a(3 

44 


A «0 


45 


A*P 
^46 


Ag Ag Ag Ag Ag Ag 


M2 


1 53 


54 


Ag Ag Ag Ag Ag Ag j 


(204) 


and the corresponding column vectors as 


{Ae*' 1 } = 


r Ae;” 1 




b i 1 

be? 




bl 

Afj’ 

Aef 

► 

and 

M - < 

% 

Aef 




b\ t 

l to? J 




bl J 


(205) 


This system can be solved by standard Gaussian elimination. However, if M subvolumes 
are included in the fiber, these equations represent a 6M system of equations, whose solution 
may pose storage problems on the computer. 
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8.2 Solution of Integral Equations by Iteration 

An alternative is to use equations 188 and 189 in an iterative fashion. As a first guess the 
integral terms in 188 and 189 can be neglected and we obtain 

D™ kl Aet? = DT jU AcJ, - 6Df jU [Ae“, - Acj] for a = 1 to M 

corresponding to the subvolumes in the fiber, and 


= -D”„Ac2, for a = M + 1 to N 

corresponding to the subvolumes in the matrix. These relations can then be substituted 
into the integral terms to yield an improved “Rayleigh-Born” approximation to D^ kl ^e* k f for 
« = 1 to N. This process can then be repeated until D "] kl converges to within a user 
specified tolerance. In essence we solve the equations 

D? jU {Aei,“} A+1 = Dfa Ac£ - 6D& [A«4 - Ac£] - 

M ±oc 

- £ (0 DZ„ r J s Q‘ «) Q" (-«) { Ae;f } A - 

/9=1 n p = 0 

N ioo 

- £ 6D° jU £££ W (0 DZ„r,fQ° «) <?"(-{) Ac?, (206) 

0=M+ 1 n v = 0 


or 


D"] u (A4“} a+i = - 6D? it , [A4 - Aefi] - 

- £ £ ( y jjfm r) J f ju klmn (r - r') rfV(r') ] {A e ;?} A - 

/3=1 Q=l\ Va y o J v 0q j 

£ «A“«£ (vjll dV(r) lll UUmnil ~ T,) dV '( r 'q 

0— M “Hi ” 


(207) 


9=1 \ ,a w: 


until the (A + l) th iterate differs insignificantly from the A th iterate. 

In the solution of the composite problem, two constituent phases, namely the fiber and 
matrix phases, have been considered. For composites with a third chemically degrading phase 
separating the fiber from the matrix, the preceding solutions may be modified by assuming 
that, in the summations from (3 = 1 to M, some of the subvolumes, say from (3 = L to 
A/, pertain to the degraded material. It will then be necessary to postulate a viscoplastic 
constitutive relation for this chemically degrading phase. 

The total strain increment in the a th subvolume in the unit cell is then obtained by 
averaging equations 57 and 80 over the a th subvolume to give, 

M ±oo 

A elr = A 4 + £ £ £ £ (0 (0 A£,A £ ;f + 

0=1 n p =0 

N ±00 

+ £ £ £ £ 9Umn (0 fQ? (4) Q 0 (-{) ^Cnr.Acf, (208) 

0=M+1 Tip — 0 
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or 


Aelr = ^l + '£ f:\^J[fdV(T)ffju u (r - r') rfV(r’) \dz„m:h + 

0=\ q=\ \ VaJ v; V„, ) 

+ tt\vJh M III a ^ T -^ dv{ r') )D™ n „Acf. (209) 

/3=M+l<z=l\ “vi Vfl, / 


for a = 1 to Ah 

The constitutive relation, required to update the stress and state variables in each sub- 
volume, is then given for the a th subvolume as the average of equation 35 in the form, 

A<r“ = DjJ„ (A4“ - Aei") (210) 

If we assume that N = 2, with one subvolume in the fiber and the other in the matrix, 
then the theory is similar to the self-consistent model in which the strain increments in 
the constituent phases are assumed to be spatially constant and equal to their respective 
constituent volume averages. However, the interaction effects of the nearest neighboring 
cells are fully accounted for since geometric periodicity is assumed in the integral equation 
formulation and the material outside the unit cell has not been smeared into an “effective 
uniform material. 

In both the Fourier series and the Green’s function formulations integrals of the form 


JJI^ KT dV( r) (211) 

need to be evaluated over the subvolume, V a . These Laue interference integrals [39] can be 
evaluated exactly if each subvolume consists of a circular or oblong cylinder. In the case of a 
circular cylindrical fiber, each subvolume within the fiber would consist of an infinite cylinder 
with a cross-section in the shape of an element of area in cylindrical coordinates, comprised 
of two circular arcs with constant radii, T\ and r 2 , and two radial segments along the lines 
of constant 9\ and 62 - An attempt will be made to evaluate equation 211 for this type of 
cross-section. If this proves too unwieldy, the subvolumes within the cylindrical fiber can 
be taken to be cylinders themselves, with the cylindrical fiber represented as a “bundle of 
sticks”. We assume that the actual fiber is comprised of subvolumes of the correct shape, but 
we make an approximation in performing the volume integration over a circular cylindrical 
subvolume. 


9 Concluding Remarks 

This document is the first annual report on NASA Grant NAG3-882. Much of the work on 
which this report is based exists only as a melange in the literature and we have therefore 
attempted to write the report in enough mathematical detail that it can be worked through 
without reference to the literature. In the second year we shall work out the required integrals 
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in the formulations and program the methods in FORTRAN subroutines suitable for inclu- 
sion in nonlinear finite element programs. In the third year we will determine the material 
constants for various composite materials and provide a comparison of the present theory 
with finite element and experimental results. 

Our aim is to produce an end product which can be used in nonlinear finite element 
and boundary element programs for analyzing the structural behavior of composite materials 
under thermomechanical loading conditions at elevated temperature. 

The viscoplastic behavior of periodic composites is analyzed by means of implicit integral 
equations. These integral equations arise when the problem of determining the stress-strain 
variation throughout a unit periodic cell in the periodic composite is solved by a Fourier 
series or Green’s function approach. In this report we show that the Fourier series and 
Green’s function approaches are mathematically equivalent by means of the Poisson sum 
formula. By applying simplifying assumptions the integral equations can be solved in an 
approximate fashion and used in structural analysis programs to analyze the overall behavior 
of the composite. When the strain-temperature history at the “damage critical” location 
has been determined from the structural analysis, this can be used to “drive” the “exact” 
integral equations to determine the stress-strain history variation throughout a unit periodic 
cell located at the critical location. 

The unit cell in the periodic structure can be formulated to analyze fibrous, laminated 
and particulate composites. By retaining the effects due to the application of displacements 
and tractions at the surface of the composite it is also possible to analyze the behavior of thin 
walled composite sections such as are found in turbine engine combustor liners and blades. 
When this is done the integral equations which must be solved are basically those which are 
used in boundary element programs. In the constitutive subroutine which we plan to embed 
in the nonlinear finite element program to analyze the overall macroscopic behavior of the 
composite, we effectively have a boundary element equation (specialized for the case of a 
periodic composite) which we solve in an approximate fashion for the stress at the Gaussian 
integration point when the boundary displacement on the element is prescribed by the finite 
element program. 

When the effects of damage are included in the constitutive formulations it will be possible 
to embed the subroutine in an optimization program such as ADS in order to determine 
optimum composite configurations. 
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Appendix A 

Properties of the Green’s Function 

Consider a point force f k (V) acting at the point r' in an infinite medium with elasticity tensor 
Dijki- From the definition of the Green’s function the displacement at the field point r due 
to the point force f k ( r ') at r' is 


Ui( r) = Gij (r - r') f 3 (r') 
so that the infinitesimal strain at r is 


£im(r) = - 


1 ( dGjj (r - r') dG m} (r - r') 


and the associated stress is 


ji dXi 


Si (**') 


or 


^fcp(r) D kp i m ^ 


1 f dGjj{ r-rQ | dG mj (r - r ; ) 


dXr 


dxi 


Si (*0 


(A.l) 


(A.2) 


(A- 3) 


(A-4) 


Since the elasticity tensor D kpim is symmetric with respect to the indices i and m, the last 
relation can be written as 


Y-' 1 _ n dG ij ( r — r') . 

&kp\ Y ) — DkpiTn o Jj ( r ) 

OXm 

For static equilibrium, we must have 

A(r') = ~ JJ n p (r)a pk (r) dS(r) 


(A.5) 


(A.6) 


where S denotes any closed surface in the infinite medium with an outward unit normal n p ( r) 
which surrounds the point of application of the point force f k (r'). An application of Gauss’ 
divergence theorem gives 




Writing 


fk ( r> ) = fj (r') JJJ 6 kj S( r - r') dV( r) 


r) (A. 7) 


(A.8) 


then gives the equilibrium requirement that 


f^jjj 

V 


r) ~ ^ i j 
t-'kprm 


d 2 G l} (r - r') 


dx m dx p 


+ 6 kj S (r - r') l dV (r) = 0 


(A- 9) 
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(A. 10) 


Since fj (V) and V are arbitrary, then 


D 


kpim ' 


d 2 G X3 (r - r') 
dXmdx p 


+ 6 k jS (r - r') = 0 


is the differential relation satisfied by the Green’s tensor function. When multiplied by fj, 
this is just Navier’s equation of elasticity with the displacement Wj(r) = G i3 (r — r') fj (r') and 
the body force set equal to <5 kj fj (r') 6 (r — r'). 

Rearranging the indices, this differential relation can be expressed as 


D ' iU ^d^~ + M(r) _ 0 


(A. II) 


The solution to the differential equation can be found by applying Fourier integral tech- 
niques. On multiplying the differential relation by e lKqXq dxi dx 2 dx 3 , i.e., by e lK ’ T dV( r), and 
integrating over all space, we obtain 


D 


*IJJ 


d 2 G kp ( r) 


JK q x q 


DC 

dx i dx 2 dx 3 + 6 ip J Jjs (xi) 6 ( x 2 ) 6 (x 3 ) e iKqXq dx x dx 2 dx 3 = 0 


dxidxj 

(A.12) 

Prom the sifting properties of the Dirac delta function the last integral is unity, so that 


r n d (aa ip ( r)\ 8 /«?*( r)\ 

+ Di>k2 d. x, 


OO 

III 

-oo 

+ D t j k3 -^ } eiKqXq dXl dx ‘ 1 dx3 + Sip = 0 (A. 13) 


Integration by parts severally with respect to xi,x 2 ,x 3 then gives 


OO f 

// 


D, jkl e lK ' x ' dG d k ^ T) 


e i(K 2 X2+K 3 x 3 ) d X2 d Xz _|_ 


Xl— — oo 


Dij k2 e 


,iK 2 x 2 

dxj 


x e i(K ' ri + K * T3) dxidx, 3 + 
0G kp { r) 


Dij k 3 C 


ik 3 x 3 dGfcp(r) 


dxj 


J X 3 = — OO 


///{ 


D/jk t 


dxj 


iKi + D 


ijk2 ■ 


d Chpjr ). 

dx 3 


J X 2 - — OC 
e i{Kixi+K 2 X2) d x \ _ 

dG kp ( i 


iK 2 + D ijk3 - Q^ iK ’^ elKqXq dxi dx ' 2 dx ' s + 


+ 6 lf) — 0 


(A. 14) 


The surface integrals are zero since dG kp (r)/dxj vanishes at the infinite lower and upper 
limits of integration, so that one integration by parts yields the result, 


-D, 


* Jlh tS 


dG kp ( r) tK , r 


e tK r dV (r) + = 0 


(A. 15) 
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p j 

A second integration by parts yields 


' ■ 

SO 

+ JJJ ^KtK.G^e^dVir) + S tp = 0 

— OO 

(A. 16) 


or 

DijkiK t KjGkp{ K) = 6 ip 

(A.17) 


where 

OO 



G tp (K) = JJ j Gk p {v)e iKr dV (r) 

— OO 

(A. 18) 


is the Fourier integral transform of G fep (r). By writing ( as a unit vector in 

the direction of 

lari 

the wave vector K, we have 



II 

^ [J 

II 

G? 

(A.19) 

- - 

in which I< = ^K q K q is the magnitude of the K vector. Then 



K 2 G tp (K) = 6 iT 

(A. 20) 

i— 

or 



K 2 D ijkl QCAp( K) = 8 ip 

(A- 21) 


The C’hristoffel stiffness tensor M (cf. [33]) is defined by the relation 



A7ifc(C) = DijklQ Cj 

(A. 22) 

— 

so that 


- - 

K 2 M ik (C)G kp (K) = 8 ip 

(A. 23) 


Premultiplying both sides by the tensor 7f“ 2 M _1 gives 


n 

GiA K) = K-'M^CQ 

(A. 24) 


The Fourier inverse of equation A. 18 gives 



oo 

Gij(r) = (2ir)~ 3 J j |e- iK - p G y (K)d 3 K 


(A. 25) 


whoie d K dl\ i dh 2 dR 3 , so that wo finally obtain the Green’s function in the Fourier 
integral form 

_1/ u r.r _ 

(A. 26) 

t n i j. \ 

~QC 

with 

Mi] \0 = (DikijQCk) 1 (A. 27) 


a (r)= [Tf d ’K r -K-r 

, JJ.I (2tt) 3 A’ 2 e 
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This representation of the Green’s function yields explicit results for isotropic and transversely 
isotropic materials (cf. Mura’s book, [24]). For cubic and general anisotropy the Fourier 
integral representation must be used. 

Often, we are concerned with volume integrals of the Green’s function and its derivatives 
with respect to r, such as Ukimn( r)- It is then advantageous to use the Fourier integral 
representation even for isotropic and transversely isotropic materials. The advantage is gained 
by reversing the order of the wave vector and volume integrations, whereby many of the 
integrations can be carried out explicitly. 

Sir William Thomson (Lord Kelvin) obtained an explicit form for the Green’s function of 
an isotropic elastic material in 1848. As an example we may deduce the Kelvin result for the 
Green’s function of an isotropic material from the Fourier integral relation. For an isotropic 
material 

Dijki = A Sij6ki + n (SikSji + SuSjk ) (A. 28) 

and so M,k(C) — ^ijkiOCj has form 

M lk { 0 = (ACiOt + fiSikOCi + KiOfc) = ((A + aOC i<* + M*) (A. 29) 


since 


«-tf + <f + *-@) + (£) + (§)-. 

The inverse tensor M ik l (C) is given by the relation 

Mi Jt 1 (0 = ((A + ^KtCfc + ^ik) — ~ ^ 

which is easily verified by showing that 


H /i(A + 2/i) 


c<* 


M~.M jk = 8 ik 


From the preceding relations 


A/,7 1 Mjk = I - 


A + /i 
H /i(A + 2/i) 


Ci C / I ( (A + n )CjCk + (i8jk 


) 


A + 1* s , c (A + £t) 2 A A (A + //)// ^ t 

' *ik ../x , 0 ..\S«S>fc _f\ , o..\GS>fc — ®ik 




fj,(X + 2/i) 


//(A -f- 2/i) 


(A. 30) 


(A.31) 


(A. 32) 


(A. 33) 


as required. 

The Green’s function is therefore obtained in the form 


G '» = ^ 3 !HYih e ' iKTiiK - < 2 ^ 3 ///; 


A + // KiKj iK 


r d 3 K (A. 34) 


fi(X + 2/i) K 4 

— OO — 

From a table of Fourier transforms (cf [36]) we find that 


- -hlllh 


e -iK.r d 3 K 


(A. 35) 
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which may be differentiated with respect to x, and x } to give 

-?/// 


= i_ f f f Ki ^- e~ iK r d 3 K 

dxidxj 


K 4 


By contracting the i and j indices we obtain 


d 2 r 1 


dxidxj 

The Green’s function may therefore be written as 

1 ( hi — 2 




Gu( r) = 


or 


G«(r) = 


where the relations 


87T)Ltr [ 

d 2 r 2 


(2^ 


dx q dx q r 


and 


(A. 36) 


(A. 37) 


2 r A + [i ^2 \ 

(AM) 

/x(A + 2/i) T dxidxj J 

„ x + fl (sij 

(AM) 

A + 2/x V r 2 J j 


d 2 r 8ij XiXj 

(A. 40) 

dx^Xj r r 3 



obtained by differentiating r - have been used. 


Appendix B 


Relationship Between Fourier Series 
and Green’s Function Approaches 

In the composite material the total strain increment Ae^(r) is periodic in r and is defined 
by the relationship 

Ae«( r) = Ae° kl + Ae w (r) (B.l) 

where Ae”, is the strain increment applied to the composite’s boundary which is equal to 
the volume average of Ae kl (r) over the unit periodic cell, and Ae*./(r) is the deviation or 
perturbation from the average value due to the presence of the fibers. 

From equations 84 and 85 the perturbed strain increment is given in the Fourier series 
and Green’s function approaches by the equivalent relations, 

1 it oo /*/*/* 

= 9««(C) ///A™.A<>V £(r - t 'W(r') (B.2) 

C n P=° J Vc 


or 

A*«(r) = fjju m (r - r') D™ r ,Ae‘, (r') dV( r') (B.3) 

V 

We now show that these equations are equivalent and that the Green’s function relation is 
the Poisson sum transformation of the Fourier series relation. 

From the definition of gkimniC) in equation 54 we may write 


or 


gu ,)(<) = 5 jCi + Vr'KXjO.) 

9 kiij( C1.C2.C3) — | (Af i 4 I (Ci>6.6)<jO + Vi 1 (Ci > Cn. C" i)CjG) 


(B.4) 

(B.5) 


where 


C = - = 


2nn t 


i/iT?)’*®’*©’ 


(no sum on i) for i — 1, 2, 3. (B.6) 


L 2 ) ' \ u 

We may therefore write 

0wi>(O =9kiij(Ci (nun 2 ,n 3 ) ,C2(ni,n 2 ,n 3 ) ,C 3 (n u n2,n 3 )) = fkiij (n u n 2 , n 3 ) (B.7) 

The perturbation strain increment can then be written in the form 

±00 ±OC ± 0 O 

Ae w (r) = -r-j—j- Y, Y Y' fw fai, n 2 , n 3 ) iff D% rs As* s (r') x 

ni=0 n2=0 n3=0 

/ 27 mj , , 27 T 712 , ,N 2jrrt3, O 

'( Lj ( x ‘ x 0 + ( X2 ^) + l 3 (* 3 (B.8) 


x e 
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or as 


where 


1 ±oo ±oo ±oo 

A £fc ,(r) = T~rr E E E'M n i,n 2 ,n 3 ) 

L\L 2 L2 ni _ 0 n2= o n3 =0 


(B.9) 


hid (ni, n 2 , n 3 ) = fmj (^1, ^2> ^3) JJJ D™ r 3 A£* ra (r') x 

/27ml, , \ 27rn 2/ , x 27rn 3 ( , x\ 

< (“Er(* l -*i) +_ Er(* a "*») + u {x3 ~ X3 >) dx\dx' 2 dx' z (B.10) 


x e 


By the Poisson sum formula (c/. Morse and Feshbach’s “Theoretical Physics , [37]) 

±oo ±oo ±00 ±oo ±oo ± L 1 L 2 L 3 f f f j 3 t^ i(miKiLi+m 2 ^ 2 ^ 2 +"i 3 ^ 3 ^ 3 ) v 

E E e hki( n u n 2, n 3) = el L w yjy dKe 

— 0 712 =0 773 =0 ^1 — 0 7712 0 7773 0 — (X) 


x h 


kl 


( 


KjLi K 2 L 2 K 3 L 3 

27r ’ 27 r ’ 27r 


) 


(B.ll) 


where the sum over the integers 711 , 712,713 is replaced by the sum over the integers mi,7n 2 ,7n 3 
in the Fourier integrals. The sum over rrii includes the case where 7ni = m 2 = m 3 = 0. 

We now have the alternative sum, 


^ ±00 ±00 ±00 

A £fci (r) = — Y, E Yl h kl {n u n 2 ,n 3 ) 

L /\ L / 2 ^s m— 0 n 2 =o 713=0 


= h "’\ 2n ' 2n ’ 2n 

TTl 1=0 771 2 = 0 771 3 — 0 _ V / 


(2tt) 3 ~ '" u \ 2 tt ’ 2tt ’ 2t r 

x JJJ D™ rs As* s (r') e K Kl ( Xl - x ') +K2 ( x *~ x ' 2 ) +K3 ( X3 - x ' 3 ')} dx\ dx 2 dx 3 (B. 12 ) 


or 


*2? *2? fff d 3 K (K\L X K 2 L 2 K 3 L 3 \ i{KlXl+K2X2+I<3X3 ) ^ 

77ll=0 77l 2 =0 7713=0 _ QO V ' 

x JJJ £)rn ^ £ ^ e -i(K 1 (x\-miL 1 )+K2(z' 2 -m2L2)+K 3 (x' 3 -m3L 3 )) ^ ^ <&' 


Due to the geometric periodicity of the unit cell we may write 

Ae* s (r') = Ae* s (x^, x' 2 , x' ) = Ae^ - imLi, x\ - m 2 L 2 , x' - m 3 L 3 ) 

and 

dxj dx 2 dx 3 = d (xj - miLi) d (x' 2 - m 2 L 2 ) d (x' 3 - tti 3 L 3 ) 
so that by making the change of variable 

(x\ — th\L\ , x' 2 — 7 n 2 L 2 ,x 3 — m z L z ) = (x", x 2 , x 3 ) = r 


(B.13) 

(B-14) 

(B.15) 

(B.16) 
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we obtain 


AEkl 


W - /// 

— no 


d 3 K 

(27r) 3 



K2L2 K 3 L 3 \ i(K lXl +K 2 X 2 +K 3 X 3 ) 

2tt ’ 2tt / 


X 



(B. IT) 


where the volume integration extends over the volume V c (mi,m 2 ,m 3 ) of the unit cell whose 
center is at the point (miLi, m 2 L 2 , m 3 L 3 ). Since mi,m 2 ,m 3 range over all integer values, 
the summation of the volume integrals extends to all the cells in the periodic lattice, i.e., 
it extends over the entire volume, V, of the composite medium. The expression for Ae^r) 
thus takes the form 


Ae w (r) 


f T f d 3 K , (KxL 1 K 2 L 2 K 3 L 3 \ tK , 

JJJ(2«r Jki n 2 tt ’ 2 tt ’ 27 t) 

— OO 

X /// D$„Ae'r, (r") <T ,K r " dV(r") 

V 


(B.18) 


By interchanging the order of the volume and wave vector integrals and noting that r" 
can be replaced by r' since it is a dummy integration variable, we obtain 

OO 

A, u (r)=///rfV(r')/// 

V -00 

Introducing (%^-, ~&r) pi 3 - 06 °f ( n i> n 2,n 3 ) in the expression for 

fkiij ( ni,n 2 ,n 3 ) = g k uj[ Ci («i, n 2 , n 3 ) , C 2 (ni,n 2 ,n 3 ) ,£3 (n x ,n 2 ,n 3 ) ) (B.20) 

then gives 


d 3 K 

(2t r) 3 


/fci 


*3 


K\Li K 2 L 2 K 3 L 3 


2t r ’ 2 tt 


2tt 


)e® D? jr ,Ae' r .(T') (B.19) 


Qklij 


Ci ( 


ATiLi A~ 2 L 2 

27T ’ 27T ’ 27T 


> C 2 


K 1 L 1 K 2 L 2 K 3 L 3 
2tt ’ 2tt ’ 2tt 


M 


KJn K 2 L 2 K 3 L 3 

2tt ’ 2tt ’ 2tt 


1 KjK k 


K 2 


K 2 


(B.21) 


with 


c* = 


Ki 


Ki 


K 

and the perturbed strain increment takes the form 


(B.22) 


Af,,(r) 




(B.23) 
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But, from Appendix A, 




which is the result obtained with the Green’s function approach. 

The Fourier series expression for the perturbation strain increment is thus identical to the 
Green’s function expression and the two are linked via the Poisson sum formula. 
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Appendix C 

Poisson Sum Formula 

In Fig. 7 the function f(x ) = [(^ - is shown on the unit cell extending from x — 0 

to x — L. The corresponding function defined on the n th unit cell to the right is given by 
f(x + nL) and the periodic function q(x), which is comprised of the functions f(x + nL ) 
defined on all of the unit cells extending from x = — oo to x = +oo, is given by 


q( x ) = Y f( x + nL ) 


(C.l) 


Each function, f(x + nL), is defined only over the corresponding n th periodic cell and is taken 
to be zero outside of the cell. Each function can therefore be represented as a Fourier integral 
and the periodic function q{x) can be written as a sum of Fourier integrals, 


q( x ) — ^2 f( x + nL) = ^2 Fourier integral of f(x + nL) 


(C.2) 


By setting x = 0 in both summations we obtain the Poisson sum formula. This method is 
outlined at the end of this Appendix. 

The Poisson sum formula can also be derived by expanding the periodic function q(x) into 
a Fourier expansion and showing that the Fourier integral sum, when x = 0, is the sum of 
the coefficients in the Fourier series expansion. 

Since q(x) is a periodic function of period L, it may be expanded into a Fourier series in 
the form 

_ i 27 rmx 

q( x ) = Y a ™ e L (C.3) 


where 


1 r L %2itmx 

a m = t e L Q ( x ') dx' 

L J o 

The object is to show that the Fourier series 


1 _ i2irmx rL j 2 k mx' 

q( x ) = 7 E e L e L Q ( x ') dx ' 


(C.4) 


(C.5) 


represents a sum of Fourier integrals. This is easily accomplished by introducing the expres- 
sion 

OO 

q (x') = ^2 f ( x> + n L) 

n— — co 

into the Fourier expansion and changing the integration variable by means of the relation 

y = x 1 + nL (C.6) 

In this way we obtain 

1 .22. j ’iwmx JE. 


. 22 , 1 ~ i2wmx °° n f 

( i( x ) = Y f( x + nL ) = t Y e 1 Y / 

n = - oo ^ m—~ oo n— -oo , . > 


i2nm{y—nL) 

e L f(y)dy (C.7) 


y=(n-l)L 
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m 


E3 




The exponential function exp(i27rm?//L) is a periodic function with period L , so that 


i2irm(y—nL) i2Trmy 

e l = e L 


(C.8) 


If we set x = 0 and note that the sum over the integration limits is equivalent to summing 
over the entire axis of x from x = — oo to x = +oo, we obtain 


y^> y^ 1 roc i 2ixmy 

?(o)= E f( nL )= E 7 / e L ~f(y) d y 

n=-oc m=— oo ^ 00 


Putting L — 1 gives 


°° °° roc 

9(0)= E /(")“ E / e ' 2lrrny f(y) dy 

^ J —oo 


(C.9) 


(C.10) 


and by changing the integration variable to K = 2iry/L, we obtain 

°° L r°o 


E /(») = E ^ l eimKL f [^UK 


(Sr) 


(C.ll) 


In three dimensions this result takes the form 

±oo ±oo ±oo ±oo ±oo ±oo T T T rYr 

E E E /(»■.»>."«) = E E E ^YlII i3K ^ m ' K ' L ' 

m=0 n 2 =0 n 3 =0 rrn- A A A 1 Z7r 1 

X 


+ 1712^2^2+^3^3^3) 


^ ^ (2rr\Z J 

=0m 2 =0m 3 =0 V Z7 V 

tfiLi K2L2 K 3 L 3 

2tt ’ 2tt ’ 2tt 


x 


'0 


) 


(C.12) 


which is the form used in Appendix B. The cubic function defined here for illustration purposes 
has the property that the constant a 0 in the Fourier expansion is zero, since f 0 L f (x) dx = 0. 
This term may therefore be omitted from the summation on the left and the summation signs 
primed to denote the omission of the term with Tij = ri 2 = n 3 = 0. 

It is now possible to show that the Poisson sum formula follows from the Fourier integral 
sum in equation C.2. 

We have the Fourier integral sum representation 


oo oo -oo 

E /( m )= E I f(y)Hy~ m ) d y 

m=— oo m—-o o 00 

where the Dirac delta function is given by the Fourier integral 

6 {y-m) = E r e i{y ~ m)z dz 

Z 7 T J —00 


(C. 13) 


Then the Fourier integral of f(m) is 


lr°° r 

/M = -/ e~ imz dz j e 

Z 7 T J - 00 J — 00 




f(y ) dy 


(C.14) 


(C-15) 
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Putting z = 2-nct gives 


/ oo roc 

e~ l2nma da I e i2iray f(y) dy (C.16) 

-OO J — OO 

The Fourier series expansion which culminates in equation C.10 shows that 

/ OO 

e' 2m *f(y)dy (C.17) 

-oo 

so that equation C.16 becomes 

/ OO 

e~ i2nma f(a)da (C.18) 

-OO 

We may therefore write 

OC — OO OO OO 

E /("*)= E/M= E /(-m)= E /_ e' 2 ""- /(a) rfa (C.19) 

m=—o c m=oo m=—oc m=— oo 00 

This is the Poisson sum formula in equation C.10, from which equations C.ll and C.12 follow. 
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Appendix D 


Integral Equation for Displacement Increment 
In Neighborhood of Free Surface 


The static equilibrium equation 
equation 37 as 

D m 

u ijki 


for a medium with elasticity tensor 



D m 

u ijki 


is obtained from 


(D.l) 


in which Ac^(r) = Ae°,(r) + Aefct(r). n , . , 

In this equilibrium relation we will not assume that the strain increment app le 

at the surface of the composite is constant and will take it as a spatial variable. If we set 


A/i(r) = -L, { D «« Ae “ (r> } 


(D-2) 


and note that 



the equilibrium equation may be written in the form 


d 2 (A uf) 

D™ \ LZ = A f. 

D « kl d Xj dx k jl 


(D.3) 


(D.4) 


where the symmetry of D™ kl with respect to the indices k and l has been used. On denoting 
the operator T it by the relationship 


the equilibrium equation is 
Now consider the integral 


Tii = m 


d 2 


ijkl dXjdx k 


TuAuf — A/j 


(D.5) 

(D.6) 


I(<t>,rP) = (D.7) 

for any two field variables 0 t (r) and 0i(r). These field variables may be tensors of any rank. 
For example, if <p and 'ip were second rank tensors, then 1(0, ip) would be a second rank 
tensor integral 

/„(<»,,/,) = ///M v')^ jq (v')dV(T’) (D.8) 

V 
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From the definition of the operator Tij we have 


«*,*)= fff*.woz,£ d 9 W-) dV{r ' ) 


(D.9) 


and since 

d 


d (rh <r’) n m = t ^ + a ( r ')D Tn — ( m 

dx'V i[ - ) Dipqj dx' } dxL ° ipqj dx' + <j) ^> tp9j dx'{ dx' ) ( 


' ™ dx' q 
the integral becomes 


10 ) 


= ///4 (D11) 

The first integral can be transformed into a surface integral via Gauss’ divergence theorem, 
so that 


W) = (D-12) 

5 <1 v P 1 

By interchanging the arguments (p and ip it is evident that 

W. *) = // «p ( r ) A M <*S(r') - JJj^p D ^,^P dV ( r ''l < D - 13 > 

5 <1 v p 9 

Now the elasticity tensor 7)™^ is symmetric with respect to its indices, so that the interchanges 

ip -»■ pi, qj -> jq, ip <-► (gj or jg), gj <-► (ip or pi) (D.14) 

leave the elasticity tensor unaltered. This shows that the volume integrals in I{<p,ip) and 
I(ip,(p) are identical, so that Green’s identity ([38], page 434) can be written as 


jj] - AKAi) dV = // n„ <<5 (D-16) 

Now choose the field variables (p and ip as a vector and second rank tensor in the forms 

<Pi (r') = A uj (r 1 ) - A«° (r') = An, (r') (D.17) 

and 

ip i:i (r') = Gij (r - r') (D.18) 

where (r') is the perturbation displacement increment Au, (r') and (r - r') is the Green’s 
tensor function satisfying the differential equation (cf. Appendix A), 


T t] G jk (r-r') + ^<5(r-r') = 0 


(D-19) 
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The integral relation then becomes 


/// {(Auf (r') - Au° (r')) (r - r') - 

V 

- G* (r - r') (Auf (r') - Au° (r')) } dV( t') 

= jj n v (r') | (Auf (r') - Au» (r')) ^ 

1 "(MW) 


- G* (r - r') D\ 


dS{v') 


(D.20) 


Replacing TijGjk (r — r') by -8 lk 8 (r - r') in the first term of the volume integral gives 
(r) - Au°(r) + JJJ G ik (r - r') T i3 (A uj (r') - A u° (r')) dV( r') 

-t! (r'Jr It-t 'ID" 

- JJ II,, (r ) jC, t (r r D ipqJ ^ ftjj j 


(Auf (✓) - Au? (r')) iS(r 


(D.21) 


From the definition of the operator we obtain the relation 

■FyAuJ(r) = A/i(r) = {D" M A£;,(r)} 

and by inserting this result into the integral equation and noting that 


(D.22) 


a(Au°(r')) _ 1 / a(Au; ( r-)) «(XW)\ _ A£0 23 

dx > q - 2 1 dx' q + dx'j I J K } 


we find that the total displacement increment is 
* T/ \ a n / \ f f f ^ t /_ _/\ ^ 


Anl(r) = Au°(r) - JJJg* (r - r') {D- ra (A^ (r') - Ae? a (r')) } dV( r') + 

+ JJn r (r') j C, t (r - r') D" , (Ae£ (r') - A^ (r')) - 
s ' 

- (Auf (r') - Auf (rO) r ' ) } iS( r>) ( 


(D.24) 


In the volume integral an integration by parts with Gauss’ divergence theorem using the 
relationship 

dGik (r — r') = dG* (r - r') 
dx\ dxj 
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gives 

jffa * (r - r') A {A", (A< s (O - Ae?, (r'))} dV(r') 

V' J 

= /// A ( G, ‘ (r - r,) D?i " ( Ae ” (r,) - Ae " (r '>)} dv(l,) “ 

V 3 

- Jff 90 *^- 0 * r. ( A < M - A 4 M) dK(r') 

V J 

= JJ n, (r') G ik (r - r') D£ rs (AeJI, (r') - Ae° rs (r')) dS(r') + 

S 

+ jjj - G ' k Q x ~ ~ D ^‘ ( A£ *> < r ') “ Ae “- ( r '>) (D ' 26) 

‘ V 3 


This result may now be substituted into the integral equation to produce 


Aul(r) - + 

V 3 

+ jj n, 00 l Gik (r - r') D", (A£ (r') - A<, (r')) + 

+ (Auf (r') - A«° M) A7„ gG,> i'~ r ' ) } « 


(D.27) 


In the first two terms of the surface integral we observe from equation 35 that 

nj (r') D? jrt (A# (r') - A< a (r')) = n, (r') = At, (r') (D.28) 

represents the incremental surface traction on the surface of the composite. Equation D.27 
represents the well known Somigliana identity ([38], page 93) for the displacement increment. 
In the case where the composite is assumed to be of infinite extent the surface integrals in 
the preceding integral equation vanish, and if Ae° s (r') is assumed to be spatially constant, 
the total displacement increment is given by the relationship, 


Auftr) = Aujfr) - Jjj dG “‘g x ^ D^.As‘, 00 WW) ( D - 29 > 

V 3 


which corresponds to equation 79 and is the form used in the main report. However, equa- 
tion D.27 must be used when the surface is not infinitely removed and if As° s (r') is not 
assumed to be constant. 

In a finite element context it will normally be assumed that the fibers are very small in 
comparison with the dimensions of the finite element. At the Gaussian integration point in 
the finite element it is then permissible to neglect the contribution from the surface integrals 
since the surface of the finite element is assumed to be many periodic cells away at “infinity” . 
In some situations, however, this may not be a valid assumption. Some turbine blades and 
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turbine engine combustor liners are fabricated from thin sections in which the central passages 
are hollow to allow cooling air to pass through the component. In the thin cross sections of 
such components the surface integrals must be retained in the constitutive formulation. 

Suppose, for example, that the total displacement increment at the node points of a finite 
element are given. From these nodal values and a knowledge of the element s displacement 
interpolation functions it is then possible to compute the total displacement increment Au£(r) 
on the surface of the element and the total strain increment Ae° 3 (r) at any point. Since 
A uf (r') = A (r') on the surface of the finite element, the last term in the integral equation 
vanishes and the total displacement increment is determined from 


AuJ(r) - Z// 90 ^ — D™, (A e‘, (r') - AeJ, (r')) dV(r') + 

V j 

+ // n, (r') G ik (r - r') A?,, (a£ M - A e ;„ (r')) dS(r') (D.30) 

in which the terms in the surface integral represent the contribution to the total displacement 
increment due to the incremental traction, 

A (. (r') = n , (r') (AeJT, (r') - Ae‘„ (r')) (D.31) 


Aul(r) = 


on the surface of the element. This surface traction is needed to maintain the displacement 
increment equality A uj (v 1 ) = A u® (r ; ), which is imposed at the element’s surface. 

By differentiating AtqT(r) with respect to Xk and xi and taking half the sum, the total 
strain increment is subject to the integral equation 


AeS(r) = Aey(r) + JJju tUi (r - r') A", (A £ ;, M - Ae°„ (r')) dV( r') + 

V 

+ II n ‘ (r '> \ fe,~ - + 9Ga £ t — ) D Zr. K- M - Ae;, M) dS( r') 


(D.32) 


in which 

M = DT, U Ac„ (r') - SD m (r') [AeJ (r') - Ac„ (r')] (D.33) 

and this integral equation should be used for thin sections of composite material where sur- 
face effects are important. This implicit integral equation is similar to that for the infinite 
medium but contains a correction term for the surface effects in the last integral. This surface 
integral will become less important — due to the derivatives of the Green’s function when 
the integration points r' are far removed from the field point r and it vanishes for an infinite 
medium. 

In the preceding development it was assumed that the displacement increment A u® (r 1 ) was 
known, by interpolation with the element displacement polynomials, from the nodal values. 
This forces the incremental surface traction Atj (r') to adopt a periodic distribution in order 
to maintain A u[ (r') = Aw° (r') on the surface of the element. We could, alternatively, assume 
that the surface traction increment is zero on the free surface of the element, in which case 


the total displacement increment Au[(r) will exhibit a periodic variation on the the surface 
and the surface takes on the appearance of a frilled structure. 

If we therefore assume that the finite element is thin (see Fig. 8); that the surfaces are 
free of surface traction; and that the surfaces at the ends of the finite element are sufficiently 
far removed from the Gaussian integration point, the first term in the surface integral in the 
integral equation is zero and in lieu of equation D.30 the relationship for the total displacement 
increment now takes the form, 

Anftr) = Au°(r) - fff ^ffe ~ (^. M - Ae?, (r')) dV( r') + 

V j 

+ JJ n, (r') (Auf (r') - Aa? (r')) — ~ dS( r') (D.34) 

s 

The solution to this integral equation gives a periodic total displacement increment, Au[(r), 
which, on the surface of the composite, will exhibit frilling. 

It is clear that during the finite element analysis frilling will not occur in the element. The 
interpolation functions normally used in isoparametric elements are linear and quadratic, and 
cannot adopt the required periodic behavior. However, the stiffness of the finite element — as 
computed at the Gaussian integration points with the composite constitutive model — will 
reflect that the fact that the constitutive properties are computed as though the element were 
free to take on a frilled appearance. When the “damage critical” strain-temperature history is 
used to determine the stress-strain history variation throughout the unit periodic cell outside 
of the finite element program, the preceding integral equation will allow the frilled appearance 
of the composite to be calculated. 
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Appendix E 


Evaluation of the Eshelby Tensor 

The Eshelby tensor S ip i m is defined by the relation 


Sinlm. ~ 


d 2 


Ufa (r - r,) dV(t ' ] + dZdx~ k JIJ Gii (r “ ^ iV{T ' ] } DiUm (E1) 


hpim 2 I dxidx k 


or as 


Sipim — JJJ Uiprs ( r — r ) dV ( r ) Drslm 

V 


(E.2) 


where the field point r lies within the volume, V, and where the volume extends over an 
infinite cylinder of radius a in a medium with elasticity tensor D ijk i. Although the Green’s 
function for transversely isotropic materials is known [24], it is more convenient to work with 
the Fourier integral representation of the Green’s function as given in Appendix A. 

Introduction of the Fourier integral representation, 


OO 

Gik (r - r') = JJJ 


d 3 K M^(C) — tK.(r-r') 
(27 r) 3 K 2 


(E.3) 


where Q = Ki/K = K t / jKqKq, into one of the volume integrals in the definition of S lvhn 
gives, on reversing the order of the volume and wave vector integrations, 


or 


Lkgij 


Lk9ii ~ dx k dx g III (2tt) 3 k 2 e 


d 2 


dx k dx g 


fffal r-r')dV(r') 


(E.4) 


' d 3 K My 3 (C) — iK.i 


"/// eKr dV ^ ( E - 5 ) 

The Laue interference integral [39] extends over the cylindrical volume and can be written as 

/ = JJJe iK r 'dV( r') = JJJ dx\ dx' 2 dx 3 (E.6) 


v v 

Let x\ = Q cos 9, x' 2 = £>sin 6. Then in cylindrical coordinates 

foo ra r2n 


i = r r r e * ^ ^ e deM 

J — oo Jo Jo 


Since 


/ OO . 

e iK * x * dx' 3 = 2nS (K 3 ) 

-OO 

where 6 (K 3 ) is the Dirac delta function, the integral takes the form 


•a /*27 t 




ji q cos 0+ K 2 Q sin 0) 


gdQ dd 


(E.7) 

(E.8) 

(E.9) 
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Let q = 


, dg = 


dq 


yfi^TW y/K? + K$ 


Then, if K = Jr? + K%, 


( qK\ cos 9 qK2sinO \ 

y/Ki+KZ + VKi+K3J qdqd 0 

y/K? + Kl Jo Jo 


(E.10) 


K K 

If we now set — ===== = cos #', ======= = = sin#', then 

s/Kf+iq Jki + ki 

j = 2 tt ($ (^ 3 ) f aK f 2 * e i qc0 s(e-6') 

Jo Jo 


\[i< I 2 + xf 


qdqdO 


(E.U) 


Since the integration extends over a whole circumference, it is immaterial where the origin 
of 0 is placed. The integral may therefore be written as 


or as 


/ = - t===== qdq [ 2 * e iqcos 

s[k\ + K\ J o Jo 

2ir6(K 3 ) J n , . 

= 7 ■■■ / qdq2nJo(q) 

sjK\ + K 2 Jo 

J 1 (a^Kl + Kfj 


dd 


I = 4ir 2 6 (K 3 ) a . 

JkiTk! 

where J 0 and J\ denote the usual Bessel functions of order zero and one. 
The integral L kg ij can therefore be written as 


(E.12) 


(E.13) 


Lkgij 


J J J (2ir 


d 3 K Mjj\ C) d 2 e~* Kr a _2 
(2i r) 3 K 2 dxkdxg 


4‘n 2 8 (K 3 ) a 


Ji (a^K 2 + Kfj 
sjK 2 + K 2 


Now 


so that 


d 2 e 


tK.r 


dxkdx g 


= -K k K g e 


tK.r 


(E.14) 


(E.15) 


r _ 1 fff dK,iK 2 dK 3 r , 

Lkg,J ~ 2n! II K\ + K\ 4- Kl KkK ’ M,i (C) 

— OO 

Jx {asjK 2 + K 2 } 


x e'^ 3X3 e ~ i{KlXi+K2X2) 6 ( K 3 ) a 


v/A"? + A'! 


(E.16) 


If k = 3 or g = 3, the Dirac delta function <5 (K 3 ) gives zero values for the integral. Hence, 
the non-zero values of L kgi j are given by k = 1,2 and g = 1,2. 


1 


m 

1 
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Invoking the sifting properties of the Dirac delta function, viz., 



f{K u K 2 ,K 3 )6{K 3 ) dK 3 = f(K u K 2 ,0) 


(E.17) 


then gives 


J, (ay/Kl + kI 


1 fisrjr »/-i ,r , r 0\( KkK ° \ -U- 

L ^ = -^jj dK ' dK ' M '> (<1.6.6 = °)i i spirKfJ e Jk f 

-CO V 


where the unit vector C, is now defined by the relations 


Ci = 


K x 


\J K\ + K\ \jKi + K% 


b. r. 


If we put 


Cl = -TT = 


K x 


= COS #, C2 = ~T7 — 


C3 = 0 


K 2 


+ AC] 
(E.18) 


(E.19) 


= sin# 


(E.20) 


(E.21) 


A^ ^AT] + AT] K yjKi + K 2 

and set x 1 = 7 ’ cos 0, X 2 = r sin 0. then in cylindrical coordinates, 

1 roc r2% . /- \ /- r 0 -iKrcos{Q-<t > ) Q^l( a AT) 

Lkgij = y o / o K dK d9 Mij (Ci,( 2 )(k( g e ^ 

The integration with respect to 9 extends over a complete circumference, so that 

L klii = -if 'm~' (C„C 2 )« a <W J™ ae~ iKrcos9 Ji(aK) dK (E.22) 

Since S ip i m is real, the real part of the preceding integral involving the integration over K is 

N = f acos(A> cos#) J\(aK) dK (E.23) 

Jo 

Setting z = r cos#, and noting that cos(A'z) = cos(-ATz), we need be concerned only with 
positive values of 2 . Now if the field point r lies within the cylindrical volume, then 0 < 2 < a. 
But, from Gradshteyn and Ryzhik [40], 


/*oo 

N = acos(Kz) Ji(aK) dK 
Jo 

If ip — sin -1 ( 2 /a), then 

AJ . a cosip 

N= 7^^ 


a cos {sin 1 (z/a )} 


fa 2 — 


for 0 < 2 < a 


(E.24) 


COS Ip cos ip 

cosip 


a 


= 1 


(E.25) 
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Thus, 

1 r 2w 

Lkgij = —— J q (Ci , C 2 ) CfcCs d9 (E.26) 

independent of position r in the cylinder as expected From Eshelby’s result. In this integral we 

have Ci = cos#, C 2 = sin0, C 3 = 0, MjJ 1 (Ci ^ C 2 ) = (CmAm jn Cn) an d ^ an d 9 are restricted 
to the values 1 and 2. The Eshelby tensor may now be written as 

Siplm = Dkjlm Mpj 1 (Cl, 6) CiCfc ^ + J q Mij 1 (Cl, C 2 ) CpOfc j (E.27) 

When C:i = 0 the ChristofFel stiffness tensor for a transversely isotropic material, My, and its 
inverse, A/,” 1 , (which applies to the homogenized medium of a composite with fibers arranged 
in hexagonal arrays) have the component forms 


M n 

= 

-O1111C1 + 

\ (Dim — D1122) C2 

M \2 

= 

M2, = 1 ( 

Dim + D1122) C1C2 

M 13 

= 

M31 = 0 


A/22 

= 

2 (E>hh - 

■ D1122) Ci + DimCf 

M 23 
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A/32 = 0 


M 33 

= 

D 1313 


Mn 1 


\ (-Dull - 

D1122) Ci 2 + DimCf 


2^1111 

(Dim — D1122) 

to i 


Af~ 1 — __ 

Dim + D1122 


-Dim (-Dim ~ Du 22) 

MJ 


A/3T 1 = 0 


M£ 


DnnC? + 

2 (Dim — D1122) Cl 


\D U n 

(Dim — D1122) 

A/23 1 

= 

A/32 = 0 


M33 1 


1 



D1313 



The Eshelby tensor can now be determined by integration in the form, 

5Dmi + D 1122 
8D m i 

= £1111 

__ 3^1122 ~ Dllll 

_ 8D im 

_ Pll33 

2 Dmi 


*^i 111 
■$2222 
<S'll22 

52233 


(E.28) 

(E.29) 

(E.30) 

(E.31) 

(E.32) 

(E.33) 


(E.34) 

(E.35) 

(E.36) 

(E.37) 

(E.38) 

(E.39) 


(E.40) 

(E.41) 

(E.42) 

(E.43) 
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$1133 = $2233 

$2211 = $1122 

„ ^ 3Dim — £>1122 

£>1212 — ‘-’1221 — off 

o^llll 

$2323 = $2332 = $1331 = $1313 = $1313 = 4 

The Eshelby tensor for tetragonal materials— which applies to the homogenized medium 
of a composite with a square array of fibers — is currently being worked out. 

The results for an infinite isotropic cylinder may be recovered by taking 

Dun = 2/i(l - i/)/(l - 2u) t D n22 = 2/u//(l - 2u), and D u33 = 2fiu/{l - 2u) (E.48) 

where // is the Lame shear modulus and u is Poisson’s ratio. For an infinite isotropic cylinder 
the Eshelby tensor reduces to 


$1111 

5 — 4l/ 
8(1-!') 

(E.49) 

$2222 

= $1111 

(E.50) 

$1122 

4v- 1 
8(1-1/) 

(E.51) 

$ 2233 
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2(1-!') 

(E.52) 

$1133 

= $2233 

(E.53) 

$2211 

= $1122 

(E.54) 

$1212 

* IT 

7 ' 

II 

CN 

CM 

1— H 

c 0 

II 

(E.55) 

$ 2323 

= $1313 — $1331 — $2332 — 4 

(E.56) 


The Eshelby tensor for both isotropic and transversely isotropic materials can also be 
deduced from equations 17.27, 17.30 and 17.31 of Mura’s book, [24], by setting g = 0 in his 
notation. 


(E.44) 

(E.45) 

(E.46) 

(E.47) 
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Appendix F 


Proof that U ijk i(x - y) = U ijk i{y - x) 

i 

From the definition of U,jki(x — y) we have 


U ijk i(x - y) = - 


1 fd 2 G ik (x- y) d 2 Gjk(x — y) 


2 y dxjdxi 


+ 


dxjdxi 


D 4 5G ifc (x-y) 5G ifc (x-y) 

But ^ 1 = £ so that 


dxi 


dyi 


d 2 G ik (x- y) d 2 G ik (x- y) 


dxjdxi 

The operator can therefore be written as 


dyjdyi 


U ijk i(x - y) = -- 

But G lk {x - y) = (^(y - x), so that 


1 f d 2 G ik (x - y) d 2 G jk (x - y) 


dyjdyi 


+ 


dyidyi 


') 


or 


. _ , 1 ( d 2 G lk (y - x) d 2 G jk (y - x) ' 

ljkl X Y 2 \ dyjdyi dy t dyi 


Uijki (x - y) = U ljk i{y - x) 


as required. 


(F.l) 


(F.2) 


(F.3) 


(F.4) 

(F.5) 
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Appendix G 


Differentiation of Singular Integrals 


In the text and Appendices we have taken derivatives of the volume integrals and written 
for example, 


I kq 


d frf dGjk (r- 

dx q JJJ dxj 

q v J 



d 2 G lk (r - 


dx q dXj 


r') 


^ D^.Ae', (r') dV(r') 
DSr. Ae‘ r , (r') dV(r') 


(G.l) 


If the integration volume V contains the field point r the integrand dG ik (r r ) jdxj is 
singular at the point r' = r, and the above operation in which the derivative is taken inside 
the integral must be treated with caution, as pointed out by Bui, [41] and Born and Wolf, 
[42]. We should, in fact, isolate a small spherical volume, D, about the singular point r - r 
and evaluate the integral according to Bui’s procedure, viz., 


q D J 

= / 1 - 
J J J OXqOXj 

’///£ PIP 3 ) -WW <”■ 11 

D q 

where we have used the fact that, if the spherical volume D about the point r is small enough 
the strain increment can be considered constant and taken to have the value at the center o 
the sphere, Ae* s (r). The integral may therefore be written as 

V-D 

-jin, (r') dS(r') < G - 3 ) 

s j 

The first volume integral is evaluated in the principal value sense as D — > 0. 

Rather than using the preceding operations outlined by Bui, we may treat G tj (r - r ) 
as a Fourier integral. The preceding operations are not then necessary and the derivative 
can be taken inside the integral. That is, equation G.l is valid when the Fourier Integra 
representation of the Green’s function is used. 
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To demonstrate the validity of equation G.l, consider the singular integral used by Bui. 
He considers the derivative of the integral 


. f 1 dt 

F{x) =L ~ log 


1 — X 


1 + X 


where — 1 < x < 1. Since the integral is known, its derivative is simply found as 


dF 


1 


1 


dx x — 1 x + 1 


(G.4) 


(G.5) 


Notice that the integrand is singular at the point t = x. Bui demonstrates that in order to 
take the derivative of the integral we must write it in its principal value sense, 


dt r di 

F( * )as te8 / 


dt 


t — x 

i=-l t= x+e 


X 


(G.6) 


and the derivative dF/dx must be evaluated by noting that both limits and the integrand are 
functions of x. Using Leibnitz’s rule for differentiating an integral whose limits depend on x 
gives 


dF 

dx 
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1 


it — x 
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X+tJ 


(G.7) 


(G.8) 


dx x — 1 x + 1 

To avoid the convected terms which arise from differentiating an integral whose limits 
depend on x , consider representing the integrand as a Fourier integral. We have, from Grad- 
shteyn and Rhyzik [43], the Fourier integral representation, 


t— — = ~i= F l/f isgn(A-) dK 

t — X y27T J—oo V 2 

The singular integral F(x) can then be written as 

- L dt wJZJ\ isM ~ iK ^ )iK 

- wJZfr^ (K]e ' K ’ dK iy K,dt 
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If we now differentiate with respect to x in the normal manner we obtain 
dF 1 f°° JT ^ fW . iKx r ‘~ iKt nl 

* “ 7^L dK i r ssn{K)e tK 




= wSJ K f2 is ^ K)e<Kx i eiK - e ~ iK ) 

- wXJ K fr^ (K) K (I+1) --‘ K(, -' > ) 


(G.ll) 


A comparison of this integral with equation G.9 shows that this Fourier integral has the 
inverse relation, 


dF _ _L 1_ 

dx x — 1 x + 1 


(G.12) 


which is the correct result. 

Thus, by expanding the integrand of a singular integral as a Fourier integral, reversing the 
integrals, taking the normal derivative, and inverting the resulting Fourier integral, we obtain 
the correct derivative of the singular integral. It is then clear that if the Green’s function is 
represented in Fourier integral form the procedure of Bui is not required. In fact, the Eshelby 
tensor in Appendix E is obtained by taking the derivative of the Fourier integral, and the 
correct result is obtained. 


Appendix H 

Origin of Self-Consistency 

Many researchers in the mechanics literature suggest that the self-consistent method has its 
origins in the present century. It would appear that the method is, however, very old and 
has its origins in the last century. In the Lorentz-Lorenz theory ([42], page 87 and [45]) of 
1880 the electric dipole moment p in a dielectric is related to the electric field E' by the 
constitutive relation p = aE', where a is the polarizability. The polarizability a is related to 
the refractive index n and the number of molecules per unit volume, N. If E is the mean or 
volume averaged field applied to the dielectric the actual field at any point is given by 

47tAT 

E ' = E + ^ p (H.l) 

where AnNp/3 denotes tlie perturbation or deviation from the average electric field. As shown 
on page 85 of reference [42] this value is estimated by smearing the effects of the molecules 
outside a spherical volume enclosing the point at which the field is observed. An analogous 
formula for statical fields had been derived even earlier by Clausius in 1879 and Mossotti in 
1850. 

Twersky [44] observes: 

In the biography of John William Strutt (third Baron Rayleigh) by his son Robert 
John (the fourth baron), the son quotes the father on the verse that faces the ini- 
tial contents page of the first four of Lord Rayleigh’s six volumes of Scientific 
Papers : “When I was bringing out my Scientific Papers I proposed a motto from 
the Psalms, ‘The works of the Lord are Great, sought out of all them that have 
pleasure therein’. The Secretary to the Press suggested with many apologies that 
the reader might suppose that I was the Lord.” The Secretary need not have been 
so apologetic. The second verse of Psalm 111 should have been augmented with 
the next three lines: “His work is honourable and glorious, and his righteousness 
endureth forever. He hath made his wonderful works to be remembered.” Depart- 
ing from King James’ translation, we may read in the Hebrew of the last verse 
of this psalm the most important of all the Rayleigh principles of mathematical 
physics — that the wise beginning of work in this field is to assume that the prob- 
lem had been considered by Rayleigh and to study his works: “The beginning of 
wisdom is reverence for the Lord; very good sense have all who do so.” 

Rayleigh [45] tackled the problem in his paper “On The Influence of Obstacles Arranged 
in Rectangular Order Upon the Properties of a Medium” and was probably the first person to 
define when the self-consistent method, viz. the Lorentz-Lorenz formula, could be expected 
to break down. At the end of his paper he states: 

The general conclusion as regards the optical application is that, even if we may 
neglect dispersion, we must not expect such formulae as (the Lorentz-Lorenz equa- 
tion) to be more than approximately correct in the case of dense fluid and solid 
bodies. 
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FIGURE 1. - TURBINE BLADE WITH PERIODIC MICROSTRUCTURE. 
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FIGURE 5. - PERIODIC UNIT CELL IN HEXAGONAL FIBER ARRAY SURROUNDED 
BY NEAREST NEIGHBOR CELLS. 


SMEARED OUT 
'EFFECTIVE' MEDIUM -7 


/ 



SMEARED OUT "EFFECTIVE' MEDIUM WHOSE AVERAGE CONSTITUTIVE PROP- 
ERTIES ARE THE VOLUME AVERAGE OF THE CONSTITUTIVE PROPERTIES 
IN TIC FIBER AND MATRIX, WHEN CONSTRAINED BY THE "EFFECTIVE" 

ICDIUM. 
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IS DEFINED OVER THE PERIODIC 


UNIT CELL OF LENGTH L, FROM x = 0 TO X = L. THE PERIODIC 
FUNCTION q(x) EXTENDING FROM - oo TO co IS OBTAINED BY SUMMING 
THE FUNCTIONS ON THE UNIT CELLS IN THE FORM 


oo 

q (x) = ^f(x + nL) * 
n = -oo 

FIGURE 7. - REPRESENTATION OF A PERIODIC FUNCTION AS A SUM OVER 
UNIT CELLS. 



9 |c THESE SURFACES ARE FAR ENOUGH AWAY FROM THE 
GAUSSIAN INTEGRATION POINTS TO IGNORE THE 
SURFACE TRACTION CONTRIBUTION TO THE TOTAL 
DISPLACEfOT INCREMENT Auj (D 

FIGURE 8 . - REPRESENTATION OF A FINITE ELEMENT FOR THIN COM- 
POSITE SECTIONS WHEN SURFACE IS FREE OF APPLIED TRACTION. 
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Abstract 

This work is concerned with modeling the nonlinear mechanical deformation of 
composites comprised of a periodic microstructure under small displacement condi- 
tions at elevated temperatures. The practical motivation for such work stems from the 
need to design and optimize new multiphase materials and to predict their micro- 
mechanical and bulk material behavior under in-service thermomechanical loading 
conditions. 

Two different methods, one based on a Fourier series approach and the other on a 
Green’s function approach, are used in modeling the micromechanical behavior of the 
composite material. These two methods are shown to be equivalent to each other via 
the Poisson sum formula. Although the constitutive formulations are based on a 
micromechanical approach, it should be stressed that the resulting equations arc 
.volume, waged jo prod*# <sv*wll “eff^ive^copstitutive r^lationf wjjich i relate the- 
: bulk, volume averaged, stress tnCfefrfent to the' bulk, Volutrti Averaged, stfamirrere- 
ment. As such, they are macromodels which can be used directly in nonlinear finite 
clement structural analysis programs. 



I. Introduction 

The ultimate objective of this work is to produce a computer program to 
analyze the heterogeneous stress and strain history variation at the “fatigue 
critical" locations of a composite structure operating at elevated temperatures. 
This paper describes some of the theoretical foundations for that program. A 
mesomechanics (Haritos et a/., 1988) approach is adopted which relates the 

r c ) 
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micrornechanical behavior of the heterogeneous composite to its in-service 

macroscopic behavior. . , . 

: '*X borftprehfeiisive application 'Of ihicromedhidics to mechanic delor- 

mation problems is given by Mura (1982) in his book Micromechanics of 
Defects in Solids. The composite materials in which we are interested have 
fibers which are closely packed together in periodic arrays. Pictures of metal 
itaatrixVcomjpos;ites (tungstep-fiberrTCwfc which exhibit a 

periodic microstructure can be found. in the article by Petrasek et al . (198 ). 

' Sorde Composites' afe actually' comprised of a periodic rnicrbstruc'trire whilst 

others are possessed of an essentially randomly distributed microstructure. 
When the fibers in a composite material occupy a large volume fraction of the 
material, the induced deformation in one fiber interacts with and alters the 
induced deformation in the neighboring fibers. When the fibers are densely 
packed the interaction effect becomes dominant and must be accounted for 
in the constitutive formulation. 

Nemat-Nasser et al. (1981, 1982, 1983) have exploited the mathematical 
simplicity of a periodic microstructure in order to develop elastic, plastic, and 
creep constitutive models for composite materials. The assumption of perio- 
dicity allows the heterogeneous stress, strain, and displacement fields to be 
expanded in a Fourier series, which greatly simplifies the ensuing computa- 
tions. This technique fully accounts for the interaction effects between neigh- 
boring fibers. Even when the composite is comprised of closely packed fibers 


inclusion of the interaction effe.cts cap. be as, or.m.ore,. important than inclusion 
of the random nature of the microstructure when the fibers occupy a large 
volume fraction of the composite material. 

The nonlinear constitutive behavior of composites with a periodic micro- 
structure can also be treated with a Green’s function approach as shown in 
the expositions be Gubernatis and Krumhansl (1975), Korringa (1973), Zeller 
and Dederichs (1973), and Barnett (1971, 1972). Here, the periodic heteroge- 
neous material property variation — due to the fibers is treated as a fictitious 
body force in the matrix material The Green’s function is used to evaluate the 
displacement due to a unit point force in the matrix material, and the actual 
displacement and any point in the composite can then be determined by 
summing (integrating) the effect due to a volume distribution of fictitious 
periodic body forces. 

Dvorak (1986) and Dvorak and Bahei-El-Din (1982, 1987, 1988) have also 
made great-progress in modeling the micromechanical behavior of nonlinear 
periodic composite materials and are embarked on a conbined experimental 
and theoretical effort. 

Work on the theoretical foundations behind the homogenization of micro- 
mechanical constitutive models to produce bulk macroscopic models has 


distributed at random,, the method gives accurate results <tviemai - in asset et 
Crl982)fbnthe“effective’ , ^lastidty:teasor/.When : depSely-p«5^bd,^bSf s ^9 I t n ^ 
i large volume fraction of the composite material, these interaction effects play 
, mi/* an/t mnct included in the calculations. It appears that 
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been underway in France by Devries and line (1987), Duvaut 

and' Mstrihdnler *(i 987T Dene 'and Leguillon ( 1982), Lene (1984, 1986), and- 


Sanchez-Palencia (1980, 1985). . 

Aboudi (1987) has recently developed a macroscopic formulation tor 
periodic composites based on volume averaging Bodner’s viscoplastic con- 
stitutive model (Bodnex, ?• 

general and is not restricted to any particular constitutive model. This work 
expands the heterOgerieons displcfcment throughout the constituent phases. - 
of the composite material as linear and high-order functions of the coor- 
dinates. Good agreement with experimental results was obtained by this 


A more general approach is adopted in the present work, where the dis- 
placement is not retricted to linear or quadratic variations throughout each 
constituent phase, but varies according to the “exact” theory of an infinite 

periodic composite. _ . , 

The purpose of the present paper is to outline briefly the Fourier senes an 
Green’s function formulations for the nonlinear constitutive behavior of vis- 
coplastic composites comprised of a periodic microstructure, and to show that 
the formulations are equivalent by deriving the Green’s function representa- 
tion from the Fourier series representation using the Poisson sum formula. 
Further details concerning the formulations can be found m a recent report 
(Walker et al., 1989). 


2. Theoretical Modeling Approaches 

A periodic composite material is supposedly acted upon by an imposed strain 
increment Ae? and responds in bulk with a stress increment toy. These values 
are then equated to the respective volume averaged quantities in order to 
obtain the “effective" constitutive relation for the composite material, i.e., 




1J 


Ao u (r) dVU) and 4a? - y Jjj 4tjM dV(t). (2.1) 

V. n 

In Section 3 it is shown that the volume averaged or “affective” constitutive 
relation for the composite material can be written as 

to? = Dy k ,to°, - y J J |A7„Ac H (r) - d%,(r) ^Ae^r) - Ac k( (r)J| dV( r), 

V c 


where V c is the volume of a unit periodic cell in the composite material, Ae*,(r) 
is the total strain increment at point r in the periodic cell due to the imposed 
uniform total strain increment As®, at the surface of the composite, and Ac k ,(r) 
is the strain increment at point r in the periodic cell representing the com- 
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, ... .posit* and -AcuMu is the. strain jpcremenf. ; at, pojnt, x, thq periodic. ceU. : 
representing the deviation from isothermal elastic behavior. The fourth rank 
tensor 5D ijtl (r) is defined by the relation 

5 D iJkl (T) = %)(Df w - D" kl ), (2-3) 

' - : where 5(r) = l in the fib£r and 5(r) = 0 in the matrix. vvith f^^denoting the - 

. . '. elasticity tensor. of thefiter.and.p^ths.t oft.hematrix.. : 

Tn the expression for the average or “effective* 1 constitutive relation m (2.2), 
the quantities Ae° k „ and <5*W r ) are g* vcn - The deviation strain increment 
Ac*,(r) can be obtained throughout the periodic cell as a function of position 
r by using an explicit Euler forward difference method, since the stress and 
state variables in a viscoplastic formulation will be known functions of posi- 
tion at position at the beginning of the increment. Everything is therefore 
known explicitly except the total strain increment Ae^ifr). 

Let the Fourier series approach described in Section 3 we find that the total 
strain increment is determined by solving the integral equation 

AeJi(r) = Ae kl — 77 X X X SkiijiQ 

K «,,=o 



{A7„Ac„(r') - SD.jjn [Ael(r') - Ac„(r')]} dV( r'), 


where the fourth rank tensor g klij (Q is given by 

<WQ = + C;C*W,7'(Q), ' ' * '* (2-5)' 

in which the Christoffel stiffness tensor M fJ (Q, with inverse is defined 

by the relation 

ms) = v-Q 

with C p = = ZJZ being a unit vector in the direction of the Fourier 

wave vector and t, = denoting the magnitude of the vector C In (2.4) 

the sum is taken over integer values in which 

2jin t 2im 2 _ 2nn 3 

{,_ T r- l, • 

and L lt L 2 , L 3 are the dimensions of the unit periodic cell in the x u x 2 , x 3 
directions, so that V c = L l L 2 L 3 . The values of n u n 2 , n 3 are given by 

n p = 0, ±1, ±2, ±3,. ...etc., for p = 1,2,3; (18) 

and the prime on the triple summation signs indicates that the term with 
n k = n 2 = n 3 = 0 is excluded from the sum. 

In the Green’s function approached the total strain increment Ae^fr) is 
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determined by solving a different integral equation, viz. 



A««(r) = A4 + 


'J){D«~ACrXr') 



-5D mnr ,{ r')[A4(r') - Ac„(r')} dK(r'), (2.9) 

where the" fourth rank: tensor {/^(r — r') given the fcl component of the total 
strain increment a( ..point r duejQ the mn component of a stress increment- 
applied at point V in the infinite matrix with elasticity tensor DZ*,, »- e -» 


U klm „(r ~r') = 


3 z G, w (r-Q ~| 

2L 5x,ax„ dx k dx m J 


( 2 . 10 ) 


and the volume integration in (2.9) extends over all the periodic cells in the 

composite material, i.e., over the entire compsite. 

The Green’s function tensor is defined in Barnett (1972, 1973) and Mura 
(1982) by the Fourier integral 


—Iff 

— oo 


d 3 K 

(2;t) 3 K 2 


e -*" - r ‘). 


( 2 . 11 ) 


in which the tensor £ is now defined by the relation (i — KJK with K 
/fCK denoting the magnitude of the vector K = (K t , K 2 , K y ). 

ir is showfifby a^plyte* fha^issonsiwa ....... . . • 

and (2.9) are identical, although the summation extends over the integer values 
n lt n 2 , n 3 in (2.4) and extends over the periodic cells in (2.9). 

' •• Both (2.4) and (2.9) are integral equations .for the determination of the to tal _ ; 
strain increment Ae&r), since this unknown quantity appears both on the 
left-hand sides of the equations and on the right-hand sides under the volume 

integrations. _ . . 

The “effective” constitutive relation given in (2.2) and the total strain incre- 
ment relation, given by either (2.4) or (2.9), contain the volume integration of 
the deviation strain increment Ac tl (r). In the periodic cell the deviation strain 
increment at any point r will be determined from a unified viscoplastic 
constitutive relation (Lemaitre and Chaboche, 1985) appropriate to the con- 
stituent phase in which the point r resides. If a constituent phase is included 
at the fiber-matrix interface, a constitutive relation can also be proposed for 
this chemically degrading phase, and the resulting inelastic strain increment 
determined for inclusion in the volume integrals. This may be important for 
metal matrix composites where boron, carbon, and silicon carbide react 
chemically with superalloy matrices at elevated temperatures. 

Equations (2.2), (2.4), and (2.9) form the basic incremental constitutive 
equations for determining the “effective" overall deformation behavior of a S 
composite material with a periodic microstructure. In order to update the 
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l*» caClTjC^T tb'b 7 |>repi£rsitibii |or' 1 - :> - r ^ 

“effective” constitutive relation over the next increment, the constitutive 

f A<r 0 (r) = b 0Jt ,’(r)iA£S(rj - Ac u (rj] (2-12) 

is used, where £> 0 „(r) = 0&r or >Dfc according as the point r is in th? fiber or 
matrix. The stress cr 0 (r) and state variables ^(r) can now be updated at each 
point r in preparation for computing Ac t ,(r) in the next increment. 

The derivation of the preceding equations and some methods for their 

solution are discussed in the succeeding sections of this paper. 


3. Fourier Series Approach 

The application of Fourier series to the calculation of the “effective overall 
constitutive behavior of periodic composites has been dealt with m detail by 
Nemat-Nasser et al. (1981, 1982, 1983). This work is used in this section to 
develop constitutive relationships for viscoplastic composite materials under 
small displacement conditions. 

The periodic composite is supposedly acted upon at its surface by a spatially 
linear displacement increment, Auf(r), given by 

Auf (r) = XjAefj + x ; Ao>°, (3- 1) 

where Ae° and Ao>° are the spatially uniform strain and rotation increments 

at the surface of the composite. . 

If the matrix material was homogeneous and had no Fibers embedded in it, 

v- tlw^trtdn wQ.vld.be homogeneous $nd giyeft.by w .. ,.. 




, d(Auf) d(Au°) 
1 dxj dx t 


Since this is constant, we may trivially volume average As? over the volume 
V of the homogeneous matrix material to obtain 

v 

which, by Gauss’ divergence theorem, may be written as 

As°- = i || i(n ; (r)Au?(r) + n f (r)Auj>(r)) dS( r), (3.4) 

5 

where the integral extends over the surface of the material and n,(r) denotes 
the outwardly directed unit normal vector at point r on the surface. Thus, by 
applying the displacement increment Auffr) in (3.1) over the surface of the 
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material to produce the surface strain increment given in (3.4), (3.2) and (3.3) 
show that the strain increment in the matrix material is spatially uniform. . 

If the displacement increment Auf(r) in (3.1) is applied to the actual com- 
posite material, the total displacement increment within the material, Au< (r), 
will vary in a periodic manner due to the assumed geometric periodicity of 
the composite material, so that 

AuT(r) = Auf(r) + Auj(r), (3-5) 

where Au°(r) is the displacement increment which would be induced in ihe 
homogeneous matrix jf the fil^r pha^ weje absen^.and AuifrJ is the pertur- 
bation or deviation from the homogeneous value dud to the presence of the 

fibers. . . . 

Corresponding to these displacement increments, the total strain increment 

at any point r in the composite, A£«(r), is given by the relation 

a,tm = ApP. a. (3.6) 


(i 


where 


1 ( d(Au%) a(Auf) 
2\ dx, + dx 


and 


Ae ti (r) = 


1 f d(Au t ) d(A u t ) 

2\ dx, dx k 


)• 


(3.7) 


with AeJ, representing the spatially constant total strain increment which 
would be produced on the surface and in the interior of the homogeneous 
matrix if the fibers were absent, and with Ae u (r) representing the deviation 
from Ihe nnifortn’ value' due tQ,thepre^en,qf thq , .• v. t 
increment AeJi(r) and the perturbed strain increment AeJ,(r) vary throughout 
the composite in a periodic manner. 

.We define the volume averaged stress and strain increments as <Aff y > and 

<A£ r >, respectively. The required “effective” constitutive equation for the 
composite material is then an expression relating the volume averaged stress 
and strain increments when these are equated to the respective values, boy 
and As®, applied at the surface. Fora function /( r), which varies with position, 
the volume average is defined by the relation 

(18) 

V 

Since the composite is assumed to be comprised of a periodic aggregate of 
identical unit cells, we may write 

< / > =^JJ|/ (r) ^ (r) . < 3 ' 9) 

v, 

where V c denotes the volume of the unit periodic cell. 


i'> 


»r : ( ■ s* ( 
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. . if we volume average the to tal strain iucrement in (3.6), we obtain . .... . . 

'Tv 1 *..-*.•* A, - "v * ■ V'V****‘ ; -» '‘C-s **ry < '.V » ► ;‘*v^ *■- - 

AeJ t (r) dK(r) = Ae£, + j|| Ae H (r) dK(r), (3. 10) 


<AeL>= hj 


(Ael) = A £ ° + <A £w >. • ‘ (3:11) 


But the volume averaged total strain increment is equated with the value 
applied at the surface, so that <A £ £> = A £& and 

.... <A£ tl > = o, -■ ’ : ; V: ( 3 - i2 > 

• ' which' sKoWS that the vblutn'e 'aVe'rage of the perturbation" strain inefement; 

Ae*.|(r), is equal to zero. 

If the elasticity tensor is denoted by D ou (r) and the inelastic strain tensor 
by e£ ( ( r), then the constitutive equation at any point r in the composite material 
can be written as 

<T./r) = A ;k ,(r)(£j,(r) - ef,(r) - oc kl (r)(T - T 0 )), (3.13) 

where 

a H (r )(T - T 0 ) = f T a£(r, T) dT, (3.14) 

JT 0 

is the thermal strain and a tl (r), a£(r) are the average and instantaneous 
coefficients of thermal expansion. 

The incremental form of Hocke’s law is 

Aff.,(r) = Dy ki (r)(A4(r) - Ac w (r)), . ^ (3.15) 

where Ac t ,(r) denotes the incremental strain representing ’the deviation from 
isothermal elastic conditions and is given by 

• A«ii(r)*= Ac^ffr) 4* aJJ(r)AT — A/kU r ) AA;kf( r )(^u( r )- — *+i( T ) &ki( r )(T ~ ToD*. 


in which the tensor AD ijU (r) represents the incremental change in the elasticity 
tensor due to the temperature increment AT. 

In a unified viscoplastic constitutive formulation (Lemaitre and Chaboche, 
1985) which is integrated by an explicit Euler forward difference method, the 
inelastic strain increment Ae£j(r) is a function of the current stress (at the 
beginning of the increment), <x i7 (r), and the current values of the state variables, 


q<(r). For example, if 


£$/ fijiprs* 


(3.17) 


then Aef- = f {j {c r rJ , q s ) Ar, and the inelastic strain increment is independent of 
the total strain increment AeJ,(r). This independence of the inelastic strain 
increment on the total strain increment is no longer true if an implicit integra- 
tion method (e.g., backward difference) or subincrementation method is used. 


cmami paqg m 

OF POOR QUALITY 






; Equivalence of Green’s Function and Series 

The elasticity tensor D ok ,(r) may be written as 


A ju( r ) — ^ijkl "b <5AjTtl( r l. 


(3-18) 


where 


<5A ;4 ,(r) = 9{r)(D[j kl - D™ kl ), (3.19) 

with 3(r) = 1 in the Fiber and $(r) = 0 in the matrix, the superscripts f and m 
referring to the elasticity tensor of the fiber and matrix, respectively. The 

constitutive equation at any point r can then be written, from (3-15), as 

. . .... . . • . Affy(r) •? (A™«.+.<5A;ki(r))<Ae? t . + 4%(d - 4<4iM),. ... (3,20). 

or 

A(J 0 (r) = D&(A4 + A £u (r)) 

- { A?« Ac k ,(r)) - 5D, jkl (r) [Ae? + Ae H (r) - Ac kJ (r)] }. (3.2 1) 

If the quantity in braces is set equal to D$,A£*,(r), that is, if 

A£iA4(r) = D;%,Ac kl (r) - SD iJkl ( r) [Ae? ( + Ae k ,(r) - Ac kt (r)], (3.22) 
then (3.21) can be written in the form 

A<r 0 (r) = A?w(A4(r) - Ae?,(r)) = A"«(A4 + Ae H (r) - A£*,(r)). (3.23) 

From the preceding equation it is evident that the eigenstrain increment, 
A£?,(r), represents the incremental deviation from isothermal elastic behavior 
in the composite material when the elasticity tensor is taken to be a spatially 
constants tensor appropriate to that of the matrix phase. 

-Newton's law for continuing. Si-atic.equUibriam’ihrbughtju^ 
crement requires that 

g(A^(r)) Q 


: dxj 


Equations (3.23) and (3.24) then require that 

d{Aa.(A4 + A£ kl (r) - Ae*»)} _ A 
— u. 


dxj 


or, if A4 is constant, 


D ra 

U ijkl 


S( As k ,(r)) 

dx, 


— n m 

— U (jkl 


^(A£?,(r)) 

dx, 


(3.24) 


(3.25) 


(3.26) 


Due to the geometrix periodicity of the composite we may expand Au k (r) 
and A£ k *,(r) in a Fourier series (Mura, 1982, Appendix 3). This gives 


(0 


±o0 ±oO ±00 

Au k (r) = X X Z' Af2 k (n,, rt 2 , n 3 ) 

n x =0 n 2 = 0 n 3 "0 


f /2jtn k 

x ex mrr 


2jin 2 

*i + — — x 2 + 


L, 




(3.27) 
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where L u L 2 , F 3 are the dimensions of a unit cell in the X lt x 2 , x 3 directions. 
The coeffipents Ad t iq the Fourier expansion are detepnined by multiplying 

each side of (3.27) by exp + ^7** + 

grating over the volume of the unit cell to give 

a «*(»!, « 2 ,n,)- 2-^- r r f j o Au ‘ (r) 

L\ L 2 A-*3 J jcj = 0 Jx 2 =0 Jx 3 = 0 

where only the terms with = n, survive in the summations. 

Equations (3.27) and (3.28) can be written in shortened form as 

Au k (r) = I f I' A*©**'. 

"f =0 

with coefTicients Aii(q) determined by the inverse relation 


(3.28) 


(3.29) 


Aii k © 


-*JJ 


Au k (r)e-‘^ r dV{ r). 


(3.30) 


where 


with 


J>. r = (x„ x 2 , x 3 X K c = LjL 2 L 3 , (3.31) 


2nrti 


?. — tL (no sum on i) for i = 1, 2, 3. (3.32) 

The strain increment A£?,(r) can also be expanded in a Fourier series to give 
• • • A£ k *,(r)-£. L T A^V 5 .', . • ( 3 ‘ 33 )- 

"p* 0 

with coefficients AeJ, determined by the inverse relation 

A£?,(E,) = -jjr J j A£*,(r)e' iVr dV( r). (3.-34) 

V' 

In (3.29) and (3.33) the prime indicates that the term with n 1 = n 2 = n 3 = 0 is 
excluded from the summations, since Au k (n t = 0, n 2 = 0, n 3 = 0) represents 
a rigid body displacement increment and A#, (n k = 0, n 2 = 0, n 3 = 0) repre- 
sents a spatially uniform strain increment. 

By substituting (3.29) into (3.7); (3.7) into the left-hand side of (3.26); and 
(3.34) into the right-hand side of (3.26), the equilibrium relationship becomes 

Bfr I I I* *(***©«, + M&itktjie*' 


±® 


= z x 

n „ = 0 


(3.35) 
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D^ J M^=-iD^ j Aem: (3-36) 

If i- JTiZ i denotes the magnitude of the vector 4, a unit vector C, in the 
direction of'^ can be written as f, = £,/£ Equation (3.36) can therefore be 
written in the form 




or 


(3.37) 

(3.38) 


eWuWAW® = -iDUA^M ' 

The second rank tensor, • v-.'. •• •?• •• •' 

AiitiQ = MJQ = KM’ (139) 

is called the Christoffel stiffness tensor and (3.38) can be written as 

eM^QAu^) = -iD^jAiW- ( 34 °) 

This equation can be inverted by premultiplying each side by the inverse 
tensor <J -2 M _1 to give the Fourier expansion coefficients 


Ad*© = ->Mr k l (QD^jA£*m~ 2 - 


(3.41) 


The expansion coefficients can now be substituted into the Fourier expansion 
of Au*(r) in (3.29) to give 


A« t (r) --SII* (3-42) 

«»®o 


.ThiSfjesuIt may hOW .^ ; sttbstit,utqd iQ(9 (3j7h.?b .that the jp^rtqrbaripn. strain, 
increment may be written as 

Ae„(r)=i i ZHc.WMjti+r 3 #; 

(3.43) 


If we define the fourth rank tensor g kUJ {Q by the relation 

g Wi (Q = 4(Mi l ©CjC. + Ma l {QM (3-44) 

then the perturbation strain increment can be written in the form 

Ae w (r) - I § I' g klii m,A£*£)e* ', (3.45) 

*P m 0 

and by inserting the relation for the Fourier expansion coefficients A££ from 
(3.34), we obtain 

A%(r) = i I t o r 9 U „K> JJJ (3.46) 

where the integration extends over the volume, V e = L t L 2 L 3 , of the unit 
periodic cell- 
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From (3,6) the total strain incremetit is given by 

I'WO J J j* D ^ Ae ** (r ' )e<V(, r ‘ <fK(r ' ) ’ (3 ; 47) . 

. which, from the definition of £>? rJ Ae*(r') in (3.22), may be written in the final 
form 

A4M - A«s + i I ? t r 9„«<0 j j j «**'" 

• . K « • • 

x {ATr,Ac„(r') - SDtjJrXAeUr') - Ac rj (r')J} dK(r'). (3.48) 

This implicit integral equation — (2.4) in Section 2 — must be solved to yield 
the total strain increment AeJ,(r) at each point r in the unit periodic cell. 

Instead of solving for Ae J kl (r) from this implicit integral equation, we could 
use (3.6) and (3.22) to eliminate AeJ,(r) from (3.48) to give an equivalent integral 
equation for A£* f (r), viz., 

D£,As?,(r) = D& Ac„(r) - 5D ijkl (r) [As?, - Ac w (r)] 

- SD iJkl ( r)i £ f Q I' 9^(Q IJj D™ „ As*(r>* dV( r'\ 

v ‘ (3.49) 

The incremental constitutive relation at any point r is given in (3.23), and 
this relation can be used to update the stress state at any point r in the unit 
pell once (3.49) is solved for Ae*,(r), Alternatively, (3.48) can be solved for A4(r) 

: i.-*- ••• : > • v a ^ ! irise'fted into (3:22) Wild (3.23). the •ov^'*%fl^Ve K 'e65lstiliitSvfr relation 
for the composite material can be obtained by averaging (3.23) over the unit 
periodic cell. This gives 

<A<r 0 > = <D£,(A4 + Ae k , — As?,)> ' (3-50) 

or 

= D-J’mAe®! + D™ k t(Ac ki y — D\ jii^As?,). (3.51) 

If we define A<r° = <Acr f/ > as the volume averaged stress increment. As?, = 
<Ae?,> as the volume averaged eigenstrain increment, and note from (3.12) 
that the volume averaged perturbation strain increment is zero, i.e. (Ac ki — 0, 
then the overall “effective” constitutive relationship is 

Ao?j = D%, Ae° u - D^Aifi. < 3 - 52 ) 

or, from (3.22), 

Aafj = D™ kl Ae 0 kt - y JJJ {D&,Ac u (r) - SD^ZAeUt) - Ac u (r)]} dV{ r). 

(3.53) 


A£«(r) — Ae t , 4- „ X X 

■ r ■ ■ ■ : *V 
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then proceeds as follows: 

(1) From h knowledge 6f the stress state throughoutthe unit periodic cell at • - 

the current time, t, calculate the inelastic strain increment A q„ r ) 

from an appropriate unified viscoplastic constitutive relation. The visco- 
plastic constitutive relation will vary according as r is in the fiber or 

matrix phase, respectively. . .. .. - 

(2) Compute the eigenstrain As£(r) throughout the umt penodic cell from 

either the implicit integral (3.49) or from (3.48) and (3.22). 

(3) Compute the stress increment throughout the unit periodic cell from (3.23) 
and update the stress; strain,- 4 nd viscoplastic state variables according to. 


the relations 

a 0 (r, £ + At) = a„(r, t) + Aa^r), 


ej(r, f + At) = ej(r, t) + A<T(r), 
q f (r, t + At) = <Ji(r, t) + Aq^r). 


(4) Calculate the overall “effective” stress and strain increment for the com- 
posite from (3.53) and update the overall “effective” stress and strain from 
the relations 

< 7 °(r, t + At) = <i”(r, t) + A<r 0 (r), 

£ ?.(r, t + At) = £°(r, t) + AeJ(r). 

(5) Repeat the preceding calculations for each incremental load step. 

The preceding algorithm makes use of the fact that the inelastic strain 
increment -A4(r) isMnde.pepde.ot.pf, the to.tal.strain - 

explicit Euler forward difference method is used to integrate the unified 
viscoplastic relations for the fiber and matrix phases. If an implicit method— 
such as backward difference or subincrementation— is used, the inelastic 
strain increment depends on the total strain increment. In this case the total 
strain increment must be obtained by iterating (3.48) in the form 


AeL(r) = A4 + yl I o j]j ** <r ~ n ^Ac rJ (r\ A£ " (0) 

v. 

-SDa T'Ktolir') - Ac„(r', A£(r'))]} JV( r'). (3.54) 

The first iterative guess can be taken as A£«(r) = Af*,, and the right-hand side 
evaluated to give an improved guess for AsJ^r). This process is then continued 

with 

(4*0}.,., - A* + £ WO {}} **' A.J.MW 

-HWOKA^O-)}, - 4c„(r’, (4«yO) Jl) <"'(0. (3.5S) 

until the Ath and (A + l)th iterates of AeJi(r) converge. ^ 


1 
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Equation (3.49) is hot so convenient for iteration as (3.48) when the inelastic 
strain increment depends on the total strain increment. It is always necessary 
to know the total strain increment in order to calculate the inelastic 

strain increment Ae^(r)). But (3.22), viz., 

D™ w A£*,(r) = D&, Ac u (r, Ae^(r)) - SD Uk ,(r) [AeJ,(r) - Ac tf (r, AeJ,(r))] (3.56) 


is an implicit equation for Ae^fr) when the iterated quantity, A£?,(r), is given. 
Equation (3.48) is therefore the appropriate equation to iterate when the 
inelastic strain increment depends on the total strain increment. For further 
details, see Walker et al , (1989). 


K..*;v 


4. Green’s Function Approach 


The equation of continuing static equilibrium for the composite material 
throughout an applied strain increment is given by 


3(A^(r)) 

dx i 


+ ATiW = 0, 


(4.1) 


where A/j(r) is the incremental body force per unit volume of the composite 
material. From (3.23) and (4.1) we obtain 


D\ 


ijkl 




(4.2) 


From this equatiqn it is clear that the .divergence of the stress variation 
procfacecf ' by Asj^(f ) may be formally * ragStrded as ^ fictitiolis bbdy force * 
increment, analogous to A/j(r), which is applied to the homogeneous matrix 
material with elasticity tensor D™ kl . The theory of elasticity for homogeneous 
materials is generally concerned with the solution of the homogeneous differ- 
ential equation (4.2) — Navier’s equation — when the right-hand side is zero. 
When body forces are present, the standard method of solution is to obtain 
the displacement solution at r due to a unit body force applied at ri . This 
solution is given by the Green’s function G i; (r — r') which gives the displace- 
ment in the ith direction at r due to a unit point force applied in the jth 
direction at r\ For a distributed incremental body force A/j(r') the displace- 
ment increment at r is obtained by summing the results for the distribution 
in the form 


Auj(r) 


J J J ^ ; (r ~ r )A//(r') dK(r'). 

v 


(4.3) 


The integration extends over the whole volume, V t of the composite material 
which may be regarded as being of infinite extent. 

When Afj(r') = Owe know that the displacement solution is AuJ(r) = Au?(r), 
corresponding to an applied uniform strain increment Ae°- on the infinite 
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‘ : b<iur.dary the h^mof-nealis matrix. For ^ 

increment given by the right-hand side of (4.2), with 4//r) u - 

for thetdtsd displacement increment Auftrfcan be wnttert as ... ..... . .. ... j, 

A»f(r) - A»? (D - J j j G„(r - ri) A(DS_AC.M) W <«> 

This corresponds to (3.5), the volume integral representing the perturbed 

displacement increment Au^r) in (3.34) an * tensor D™ kt the Green’s 

For a material which is homogeneous with elasticity tensor D, kt 

function satisfies the differential relation (Mura, 1982, p. ) . 

-- -- 

by the relation /aa\ 

«5(r _ r ') = 5{x, - x\)S{x 2 - x' 2 )S(x 3 - x 3 ). (4-6) 

By applying Fourier integral transform techniques the Green s tensor is shown 
(Barnett, 1971, 1972) to have the Fourier integral form 


f d 3 K M;J l 

W ~ r ’ ] = J J J (2ir) 3 K : 


i7 l (C) 


in which the inverse Christoffel stiffness tensor MjHQ is defined by 

' MJ‘©-(W. (48) 

• , r k l IK K — K IK being a unit vector in the direction of the 

Fourier wave vector K™and K = denoting the magnitude of the wave 

vector K. 

Making use of the relation 

G ik (r - r')^;(D^ n AeL(r')) = “ r ') D «-» A£ - (r 


■ - G,t i r , - D ?, .ACn (4-9) 


we may write (4.4) in the form 


Auf(r) . A»?(.) - j{{!| (G.(r - OH-W W) 

V 


dx { 


C -3 
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The first volume integral can be transformed into a surface integral via Gauss’ , ■ 
divergence theorem, viz., 


dx\ 


7 (G iJt (r-r')Dr? mn AC(rf)dF(r') 


t/J 


n l (r')G iA (r- r ')D k n ? mJl A£: rt (rf)dS(r'). 


(4.11) 


The surface integral extends over the entire outer surface of the “infinite” 
matrix material. Since this is assumed to be at an infinite distance, all the 
integration points r' in the surface integral are at an infinite distance from the 
field point r and G a (r - r') = 0. Thus, for an infinite body the first volume 
integral in (4.10) vanishes. This would not be the case for a finite body in which 
the field point r is close to the surface integration point r', and the volume (or 
surface) integral would need to be retained for these situations. In this case 
other surface integrals would arise (Korringa, 1973; Walker et ai, (1989) due 
to the application of boundary incremental displacements or surface tractions 
on the surface of the material. 

From the properties of the Green’s function. 


dG ik (r - f) _ dG ik ( r - f) 

dx\ dx t 

wfiich follows since G ik is a function of ’ ’ ' 

r — r' = (Xj — x' lt x 2 - x’ ly x 3 - x 3 ). 
Equation (4.10) may then be written alternatively as 


(4.12) 


(4.13) 


Au7(r) = Au?(r) - 


rr«M 0: 

J dx, 


H..AC(r')^(r'). (4.14) 


But AeJ(r) = $(d(Au7(r))/dx; + 5(Au/(r))/5x i ), so that by differentiating (4.14) 
with respect to x, and Xj and taking half the sum, we obtain 


AsT(r) = As?- + J 


U ljtl (r - T')D? lmn AeUr’) dV{?), 


(4.15) 


which, by means of (3.22), may be written as 

AsJ(r) = As? + JJJ U tj Jr - f) 


-SD Un 1r')[A&f) - Ac„(r')]} dV( r'). 


(4.16) 
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An equivalent integral equation, invbl ving the eigeftstraiis incrtmeiit Ae$ (r), 
can also be obtained by using (3.22) to eliminate AsJ(r) from (4.15), which gives 


D£,AE?,(r) = D^ kl Ac k! (T)-5D iJkl (r)lAe^ - Ac tt (r)] " 


'Ill 


- SD uu (t) U k U r - ?)DZr,&eW) dV{ ft (4. 1 7) 


In the preceding equations the operator 


U ijkl (r - r') = 


l/ ^G 0t (r-r') d 2 G jk (r - r'A 
2 \ dxj dx t dx K dxi ) 


(4.18) 


gives the ij component of the strain increment at point r due to an applied 
stress increment component kl at point f in an infinite homogeneous medium 
with elasticity tensor D™ kl and Green’s function given by (4.7). 

From (3.6), (3.46), and (4.15) we see that the perturbed strain increment, 
Ae fc/ (r) = Ae kl (r) — As*,, is given by the equivalent relations 


Ae*,(r) = -^ Z Z Z' 9ktmn(Q 

K o 


III- 


,Ae*(r')e*< r - r '> dK(r'), (4.19) 


v< 


or 


As kl (r) = 


U klmn (r-r')DZrsW,(rldV(r'y 


(4.20) 


• The' volume integral in the Fourier series representation extends over the 
volume, K c , of the unit periodic cell and the summation extends over the 
integers n p = 0, ±1, ±2, . etc. y where p = 1, 2, 3. In the Green’s function 
approach the volume integral extends over the entire infinite medium, i.e., over 
all the periodic cells comprising the material. It is shown in Section 5 that the 
Fourier summation expression in (4.17) can be converted into the Green’s 
function expression by means of the Poisson sum formula. 

From (3.22) it is evident that if the elastic properties of the fiber are the same 
as that of the matrix, then 5D iJkl ( r) = 5(r)(D^ r - D™ kl ) = 0, in which case 

AfS(r) = Ac u (r) (4.21) 


is known explicitly without having to solve the integral equation. From (3.48) 
and (4.16) it can also be observed that AeJ,(r) is known explicitly when 
5D iJkl ( r) = 0. The explicit relation in (4.21) holds only when an explicit Euler 
forward difference method is used to integrate the viscoplastic constitutive 
relations. For implicit integration methods in which the inelastic strain in- 
crement Ae£ f (r) depends on the total strain increment Afi^(r), (3.48) and (4.16) 
show that even when SD ijkl (r) = 0, the equation to determine Ae^(r) is still an 
implicit integral equation. 
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5. Relationship Between Fourier Series and Green’s Function Approaches 

In the composite material the total strain increment is periodic in r and 
is defined by the relationship 

A £ L( r ) = A £ kl 4- A£ w (r), (5.1) 

where Ae% is the strain increment applied to the composite’s boundary and is 
equal to the volume average of A£^(r) over the unit periodic cell, and A£ fcf (r) 
is the deviation or perturbation from the average value due to the presence 
of the fibers. 

From (4. 19) and (4.20) the perturbed strain increment is given in the Fourier 
series and Green’s function approaches by the equivalent relations 

Ae w (r) = % 9k,ljiQ dV{r ' ] ' (5 ' 2) 


or 


K 


— in 


U kltJ (r - r’)D™„Ae*(r') dV{r’). 


(5-3) 


We now show that these equations are equivalent and that the Green’s func- 
tion relation is the Poisson sum transformation of the Fourier series relation. 
From the definition of g klm ,(Q in (3.44) we may write 

ffui/Q = i(M£ W« + Mi‘(QCA). (5-4) 

or 

g k tij(C £2* £3) = 2(^1 * ! (£i, C2J CalCjCt + 1 (C 1 » £2* £3 )£;£*)» 

■ where- ' • ..v •••/■■• ' •*’ 

2 7CFlj 


Sf — e ~ 




(no sum on i) 


Z 2 Jin, 

for i = 1, 2, 3. (5.6) 

We may therefore write 

g k tij(Q = g k tij(C l( n l> n l » ”3)* £2(^1* n 2» n 2» ”3)) = n l> n 3^> ^ 

and the perturbation strain increment can then be written in the form 

1 ±00 ±ao ±«0 f f f 

A%(r) = 7 — : — 7— E I T fkiiji n i>n 2 ,n 3 ) I Ay^Ae*^) 

LiLjL-s n,*0 I» 2 “0 i« 3“0 J J J 

V € 

K lnn. 2nn 2 

-^-i(Xi - X,) + -J— (x 2 - x 2 ) 

+ dx' 2 dx' 3 , (5.8) 
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or as 


where 


Ae u (r) = Z Z T' M"i» "i* "a* <«> 

X-n n j = 0 nj*0 n j = 0 

kkt( n 1* n 2» n 3) ” Aly(^l» n 2> ^3) ^ ^ ^ A;rs^ £ r-j( r ) 


f r27in, , .. 27rn 2/ 

x exp ji -j- L (x l - x x ) + -j—( x 2 - x 2 ) 


+ T7 (X3_ ^ ) } dx 'i dx 'i dx >- < 5 - 


10) 


By the Poisson sum formula (Morse and Feshback, 1958) we may write 

±00 ±00 ±00 

I [ r K,(n u n 2 , n } ) 

r»i =0 rt 2 = 0 n 3 = 0 


±» ±oo ±oo T 7 T 

Y V V 

L h L 

m , = 0 = 0 m 3 =0 


oo 



“X 




X /l fcl 


2n ’ 


*2±2 gjM 

2rt ’ 27t y 


(5-11) 


where the sum over the integers n lt n 2 , n 3 is replaced by the sum over the 
integers m lt m 2 , m 2 in the Fourier integrals, the sum overm,- including the case 
where m 2 = m 2 = m 2 = 0. 

We now have the alternative sum 


| ±30 ±X ±X 

Ae„(r)= - II T "i.«s) 

L>lLi 2^3 « k =0 « 2 *0 « 3 **0 


d 3 K 


» K i L t 4-m 3 K i 


x 

±00 ±00 ±ao f f 

= Z Z Z p_\3 ' 

m t = 0 m 2 = 0 m 3 =0 J J J W 

-co 

(K l L t Kil y\ 

* Juij \ 2n ’ 2jt ’ 2n J 

xJJJ P? j Ap» (r , )f» < l |: »< 3[ »~*'«>* K > (jr 7~*'» ) ' l ' l: J <3tj ~ x ' il1 dx\ dx 2 dx' 3 , 


y. 


(5.12) 
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or 


*5? t? ff fd 3 K 

A£ ‘ ,(r ) ” m^O m^O j J J (2?t) 3 

— 03 

f (K 1 £- j KiL 2 K3L3 \ „i(JC,X, +KjXj+KjX 2 ) 

* fkUi { 2n ' 2n ' 2n ) 


( 


DZ, Ae*(r') 




x e -^i(*'«-*"iL|) + ^2(*i“m 2 L 2 ) + C,(xi-m J L J )] ^ J x ^ (5.13) 

Due to the geometric periodicity of the unit cell we may write 
Afi*(r') = Ae*(Xi, x 2 , x 3 ) = Ae*(xi - m t L t , x 2 - m 2 L 2 , x 2 - w 3 L 3 X ^ 

and 

dxi dx 2 dx 3 = d(xi — m i L l )d(x 2 — m 2 L 2 )d(x 3 — m 3 L 3 ), (5.15) 

so that by making the change of variable 

(xi — m i L l , x 2 — m 2 L 2 , * 3 — ^ 3 L 3 ) = ( x i» x 2 » x 3) = r » (5.16) 

the perturbation strain increment is 

* f f f ^2^2 i(X 1 x l +X 2 x 2 +X 3 x J ) 

A£w(r) = J J J w? fklii vir • "5T* J 

“GO 

±od +00 ±oo r r r ,, 

’‘.kJU. JJJ (5 , 17) 

J%(m 

where the volume integration extends over the volume V Q (m x , m 2 , m 3 ) of the 
unit cell whose center is at the point m 2 L 2 , m 3 L 3 ). Since m x , m 2 , m 3 

range over all integer values, the summation of tthe volume integrals extends 
to all the cells in the periodic lattice, i.e., it extends over the entire volume, V, 
of the composite medium. The expression for A£ H (r) thus takes the form 

, „ rfr^K, ( k 2 l x k 2 l 2 
A£k ' r “ J J J (27i) 3Jr ^\ 2n ' 2n ' 2n ) 

D^,Ae*(T")e- X r ~ dV{r"). (5.18) 


By interchanging the order of the volume and wave vector integrals and 
noting that r" can be replaced by r' since it is a dummy integration variable. 
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we obtain 
Aa 


-IfHIfiS 

y -<n 

/ K t L t KjLj (5.19) 

A,y V 2n ' 2n ' 2n J 

Introducing [K,LJ2n, K 2 LJ2k, K 3 L 3 /2k) in place of (n l5 n 2 , n 3 ) in the 
expression for 

/ k , 0 (n l ,n 2 ,n 3 ) = ff wv (C,(ni,n 2 ,n 3 ),C2(n i .”2. n 3).C3( n i* n2 v n 3^ ( ’ ) 

then gives 

K l L l K 2 L 2 K 3 L 3 \ K 3 L 3 ^ 

2n ’ 2 k ' 2k )' 

')) “ + ^P- K J K )’ (5 21) 


Cj 


"A 1,1 \ 2k ' 2k ’ 2 w 
K ! L l K 2^2 ^3^3 


2 7i ’ 27r 2 jt 


with 


C, = -F^ 


K V ^’ 

and the perturbed strain increment takes the form 

oo 


(5.22) 


AC|u(f) — 


dV( r') 


( 2*) 3 

V - -m 

X e X v - r) D^ t Ae*{t'). 


f d 1 / ^ K + M * ®KjK k 
2 \ K 2 ' ' K 2 J 


(5.23) 


But, from (4.7), 


f f f d 3 K Mr k l ( 

G ifc ( r -n = 


(0 -iK(r-r) 


ff fd 3 K AC(Q , K(r -n 

- JJJ ( 2 k) j * 2 


(5.24) 


since G ik (r — r') = G ik ( r' — r), and therefore 

a 2 G jk (r-0 


G ik (r - r') _ f f f<f 3 K Mj?(Q K K piK . (t - n (5.25) 

JJJ ( 2") 3 K 2 ' 

— 00 

Inserting the last relation into the expression for Ae w (r) then shows that 
rrr ( 1 fd 2 G ik {r - F) 3 2 G il ( r - r , )\ „„ . / 


f 1 ft 

Ae k |(r) = - ^(r')-(- 

J J J ' 


d K(r') - 1 - T V \ — + ~ Ae* (O- (5-26) 

v 2V dx^dx, 5x y dx* J 
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From the definition of the tensor l/ u#W! (r — r') in (4.18), we see that 


A£ - (r) =JJf 


U ktlJ (T - r')D“„A£*(r') dV(r'), 


(5.27) 


which is the result obtained with the Green’s function approach. 

The Fourier series expression for the perturbation strain increment is. thus 
identical to the Greeen’s function expression and the two are linked via the 
Poisson sum formula. 


6. Concluding Remarks 

The Fourier series and Green’s function representations have been shown to 
be equivalent approaches by means of the Poisson sum formula. This method 
is well known in mathematical physics and is used extensively to turn slowly 
convergent Fourier series into a series of rapidly converging Fourier integrals. 
Both representations offer promising approaches to modeling the viscoplastic 
behavior of metal matrix composites at elevated temperatures. Having shown 
their equivalence we are free to choose between them based on mathematical 
and/or numerical convenience. Each is expected to be suited to different 
situations with respect to convergence of the series with increasing fiber 
volume fraction. Future work will explore the relative advantages of each 
formulation and the overall usefulness of these approaches in modeling the 
nonlinear viscoplastic deformation behavior of metal matrix composites. 
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Abstract— Local elastic fields in the unit cell of a periodic composite are examined numerically 
with an integral equation approach. Techniques of Fourier series and Green's functions are used 
to construct the integral equations. Numerical solutions are obtained using the Fourier series 
approach with rectangular subvolume elements. Specific results are given for a tungsten/copper 
metal matrix composite. 


1. INTRODUCTION 

The combustion chamber in the three main engines of the space shuttle has a liner material 
which is fabricated from a copper alloy. Temperature gradients are generated within this 
liner material during the space shuttle’s launch which are large enough to cause sub- 
stantial amounts of thermally-induced deformation. A tungsten fiber/copper matrix 
(W/Cu) composite is being considered as a substitute to increase the strength and improve 
the durability of the combustion liner, and may be characterized as a ductile/ductile-type 
composite material. 

Prediction of the durability of continuous-fiber-reinforced metal matrix composites 
requires an understanding of the dominant failure mechanisms in such materials. A 
requisite precursor to this understanding is the ability to predict the overall structural 
response of the combustion liner in a finite element code. Since the tungsten wires have 
diameters of about 0.2 mm, it is clear that a finite element mesh sufficiently fine to 
delineate the deformation behavior in and around the fibers on a local level is prohibitive. 

A structural analysis under thermomechanical loading conditions is feasible if the 
composite can be replaced with an equivalent homogeneous material which has the same 
overall stress-strain (constitutive) response. Armed with the homogenized constitutive 
relation, the structural analysis can be used to locate those points in the component — the 
damage-critical points — which experience the largest stress-strain excursions throughout 
the applied loading history. The strain and temperature histories at the damage-critical 
locations can then be used as boundary conditions on a small volume element to deter- 
mine the local stress, strain and temperature field histories in and around the fibers. These 
fields can then be used to estimate the durability of the component. In this paper we 
develop incremental constitutive relationships suitable for the nonlinear viscoplastic 
solution of the local stress-strain behavior. These are then specialized in numerical prob- 
lems to obtain the local elastic response in a fibrous W/Cu composite. 

2. LOCAL AND HOMOGENIZED RESPONSE 

In order to perform a structural analysis of a fibrous composite component, it is 
necessary to divide the structure into finite elements, one of which is shown in Fig. 1 . Point P 
in element ABCDEFGH represents one of the Gaussian integration points at which the 
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constitutive response is used to generate the stiffness matrix of the finite element. Ideally, 
it would be desirable to use elements that are much smaller than one of the unit cells, 
QRST, of the periodic composite,- but this would tax computer resources. Instead, if the 
volume-averaged, or homogenized, constitutive properties surrounding the Gaussian 
integration point P can be calculated, these properties can be used to compute the stiffness 
of the finite element in the structural analysis. Once the strain-temperature histories at the 
damage critical locations of the structural component are established from the finite 
element analysis, these histories can be imposed at the nodes in element ABCDEFGH and 
used to determine the local stress-strain state in the typical unit periodic cell QRST by 
means of a Fourier series or Green’s function approach (Walker et al> 1989, 1990). As far 
as the Gaussian integration point is concerned, the surface of the finite element is con- 
sidered to be many unit cells away, so that the problem of determining the local fields 
within the unit cell reduces to determining the response within a periodic cell of an infinite 
lattice when the strain increment given by the finite element code is applied at infinity. 

We therefore attack the problem in two ways. 

First, a Fourier series or Green’s function method is used to determine the 
stress-strain variation throughout the unit cell, QRST, when a known strain increment, 
say Ae° kI , is applied to the nodes of the element ABCDEFGH. This is equivalent to the 
problem of determining the local response at any point r within the unit cell of a periodic 
lattice when the total strain increment, Ae*,, is applied at infinity. The local response at 
any point r within the unit cell is obtained from the relation 

AeJ/(r) = M klrs { r) Ae° , (!) 

where M kirs { r) represents the magnification or strain concentration factor that magnifies 
the strain increment applied at the surface of the finite element i.e. at its nodes and 
gives the strain increment at any point r in the unit periodic cell, QRST. The tensor mag- 
nification factor M k]rs ( r) is a complicated function of the geometry and constitutive 
properties of the constituent materials comprising the unit periodic cell which has dif- 
ferent, but mathematically equivalent, representations in the Fourier series and Green s 
function approaches (Walker et a/., 1989, 1990). Once the total strain increment Ae t/ (r) 
at any point r is known, the stress increment can be computed via Hooke’s law in the form 

A< 7 y(r) = D ukt {v)(Ae k{ (r) - Ae^r) - a*/(0 AT(r)), (2) 

where at the point r, Djj kl ( r) is the elasticity tensor, A£*/(r) is the inelastic strain increment, 
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and «w(r) AT(r) is the thermal strain increment. The inelastic strain increment can be 
computed explicitly at the point r because the stress is known as a function of P os,t '°' 1 0 r 
at the beginning of the increment. The overall, or homogenized stress increment, A a u , 
required for calculating the stiffness of the finite element, can then be obtained by volume 
averaging over the unit cell in the form 


= - 




AOijix) dF(r), 


(3) 


where V. denotes the volume of the unit cell, QRST. , _ . 

Second, once the homogenized stress increment, Aa, jt is calculated at each Gaussian 
integration point in each finite element in the composite structure, the finite elan 
analysis will yield the strain-temperature histories at the damage cnttcal locations These 
strain-temperature histories can then be applied incrementally to the ““ 
ABCDEFGH containing the damage-critical Gauss point, and the Fourier series o 
function methods will yield the local variation of the total strain increment from ( )• 

It may thus be seen that the methods are used in a complementary fashion First to 
homogenize and obtain the overall macroscopic response of the composite, and then to 
“zoom in” and calculate the local response in and around the fibers in a unit periodic «»• 
In obtaining the overall homogenized response it is necessary to use rapid methods for 
estimating the magnification tensor M ljkl (r), because this is used at each Gauss point of the 
structure for each strain increment of the loading history. A much more accurate value of 
the magnification tensor, r), can be used in postprocessing the finite element results to 
look at the local stress-strain variations throughout the unit cell. 

2. 1 . Homogenized macroscopic equations 

It is supposed that the periodic composite material is acted upon by an imposed strain 
increment Ae? and responds in bulk with a stress increment Ao (> .These values are then 
equated to the respective volume-averaged quantities m order to obtain the effective con- 
stitutive relation for the composite material, i.e. 


*4 = v 


Aa„(r) dF(r) and 


*4 = y 


AeJ(r) dF(r), (4) 


where V is the volume of the body. . , 

The volume-averaged or effective constitutive relation for the composite material can 

be written (Walker et a!., 1989, 1990) as 

A of, = Df Jk , Ae° kl - i [ [[ \D? Jkl Ac*,(r) - dD iJkl (r)[ As T kl (r) - Ac„(r)]J dV(r), (5) 


K__ . 

where V is the volume of a unit periodic cell in the composite material, AeJ,(r) is the total 
strain increment at point r in the periodic cell due to the imposed uniform total strain 
increment Ae° kl at the surface of the composite, and Ac* ; (r) is the strain increment at point 
r in the periodic cell representing the deviation from isothermal elastic behavior, i.e. 

Ac kt (r) = Ae^(r) + a kl (r) AT(r), ( 6 ) 

where Ae^(r), a„(r) and AT(r) are the plastic strain increment, the thermal expansion 
coefficient, and the temperature increment at point r. The fourth-rank tensor 5D ljkl { r) is 
defined by the relation 

SD ijkl ( r) = d(x)(D[ kl - DT jkl ), W 

where 0(r) = 1 in the fiber and d(r) = 0 in the matrix, with D f iJkl denoting the elasticity 

tensor of the fiber and D™ kl that of the matrix. 

In the expression for the average or effective constitutive relation in (5), the quantities 
Ae° DT, kl and SD ukl ( r) are given. The deviation strain increment Ac„(r) can be obtained 
throughout the periodic cell as a function of position r by using an explicit forward- 
difference method because the stress and state variables in a viscoplastic formulation will 




32 K. P. Walker et at. 

be known functions of position at the beginning of the increment. Everything is therefore 
known explicitly except the total strain increment Ae*,(r). 


By ai 
that (8) ai 
n 2i n 2 in 


2.2. Fourier equation overview 

In the Fourier series approach we find that the total strain increment is determined 
by solving the integral equation (Walker et al ., 1989, 1990) 

Af£,(r) = A e% + — £ E L SkiijiQ 

y c n p = 0 

jj Ac„(r') - dD Urs (r')[Ael(r') - Ac„(r')]] dK(r'), (8) 

where the fourth-rank tensor gkujiQ gb/en by 

gk uA 0 = W,Mr k l (Q + CAM.7%)), ' (9) 

in which the Christoffel stiffness tensor A/y(Q, with inverse is defined by the 

relation (Barnett, 1972) 

M„(9 = d°) 

with Cp = (jfp/V^m = <?p/^ being a unit vector in the direction of the Fourier wave vector 
and £ = denoting the magnitude of the vector In (8) the sum is taken over 

integer values in which 


*i = 


2 nn x 


1 Li 1 


2nn 2 r 27 t/z 3 

<s= ir- 


(id 


and where Lj , L 2 , Z, 3 are the dimensions of the unit periodic cell in the x u x 1 , x 3 direc- 
tions, so that V c = L,L 2 L 3 . The values of n lt n 2 , n 3 are given by 

n p = 0, ±1, ±2, ±3 etc., for p= 1,2,3, (12) 

where the prime on the triple summation signs indicates that the term associated with 
- n 2 = - 0 is excluded from the sum, 

2.3. Green's equation overview 

In the Green’s function approach the total strain increment A^(r) is determined by 
solving a different integral equation (Walker et aL y 1989, 1990), viz. 


Afil/(r) = &e° kl + ~ r ') 

x [DZnrs Ac„(r') - <5D ffl „ rj (r')[A£(r') - Ac„(r')]j dF(r'), (13) 

where the fourth-rank tensor U klmn (r - r') gives the kl component of the total strain 
increment at point r due to the mn component of a stress increment applied at point r in 
the infinite matrix with elasticity tensor i.e. 

/. _ a 2 n. ( r - r'l\ 

(14) 


(W r - r') = -- 


1 fd 2 G km ( r - r') 3 2 G, m (r - r') 


3x, 


3** 3x„ 


and the volume integration in (13) extends over all the periodic cells in the composite 

material, i.e. over the entire composite. 

The Green’s function tensor is defined by the Fourier integral (Barnett, 1971, 1972, 
Mura, 1987) 


G,/(r - r') = 


oo i 3 


d 3 K M ~ *(Q r K (r - r ') 

-oo (2tt) 3 A' 2 


(15) 


in which the vector (; is now defined by the relation (, = K t /K with K - \ ! K q K Q denoting 
the magnitude of the vector K = (A) , K 2 , K } ). 
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denoting 


n 2 , n 3 in (8) and extends over the periodic cells in (13). 

2.4. Integration of the equations 

S S:S. fiber-matrix interface, a constitutive relation can also be propo ed 

OVe f q Sr"a“(.3) form the basic incremental constitutive equations for 

phases in preparation for integrating the effective constitutive relation over the 
ment, the constitutive relation 

Affy(r) = D iJkl (T)(&£ T kl (r) - Ac*,(0), (16) 

a i, n M-n f or D m t , according to whether the point r is in the fiber or 
matrix.' This reiaSn is usedm update the stress „„(r) and. in turn, the internal viscoplariic 
slate variables «,( r) at each point r in preparatton for eon.put.ng Ac„(r) at the 

'“derivation of the preceding equations and some ^V^fThe^urie' 
discussed in Walker et a!. (1989, 1990). Some numerical elastic solutions of the Jour 
series integral equation for Afil(r) are obtained in the remaining sections of the paper. 

3. NUMERICAL SOLUTION OF INTEGRAL EQUATION 

Determination of the stress and strain increments throughout the fibrous composite 
material under isothermal elastic conditions requires the solution of the integral equal 
(8), which reduces to the two-dimensional form 


Ae*,(r) = A £ °i 


. ±oo 

— ££' Skin 
Ac n„ = 0 


,(Q 


e^ u-f’) SD mnrs (r') A£(r') dS(r'), (17) 


where A c - L,L, is .he area of .he uni. ceil. a^ »^re .he .wo-dimensiona, Fourier sum 
ranges over the integer values «, = 0 , ±1, .... ±~ and n^-v, -J. • - - 
prime on the sum indicating the omission of the term m which in, " 2 i ^ 

Nemat-Nasser and his colleagues (Nemat-Nasser and Taya .1981 Nemat Nassei et , 
1982- Nemat-Nasser and Iwakuma, 1983; Iwakuma and Nemat-Nasser, 1983) have 

demonstrated .hat good accuracy can be achieved by di ; idin f lh ' 
of subvolumes, where 4s'(r') in the /i.h subvolume integral is replaced by 


AcJ,(r') - Aejf = 


AeJjfr') d5(r'). 


(18) 




which corresponds to its average value in the 0th subvolume whose cross-sectional area 


is A 0 . 
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Let there be N subvolumes in the unit cell, with M subvolumes in the fiber and 
N - M subvolumes in the matrix. Then the preceding integral equation can be written as 

AeJ,(r) = Ae° kl - f I e-'^' dS(r')^ nri A e J/, (19) 

0=1 Tip — 0 J J A(j 

where dlf mnn = - D^ nrs or according to whether the subvolume P is in the fiber 

or matrix, respectively. 

If we use Nemat-Nasser’s notation and write 


G"(0 = 


1 




dS(r), 


then the preceding equation may be written as 

Ael(r) = Ae° kl - £ £ £ ' g klmn (Q AeJ/, 


0 = I n„=Q 


where 


f 0 = 


( 20 ) 


( 21 ) 


( 22 ) 


U1 

gl 

tc 


tl 


is the volume fraction of the y?th subvolume. We may now volume-average (21) over the 
ath subvolume to obtain 

Ae T k f = Ae 0 kl - £ fS^AeJf, 

0=1 

where 

±00 

= ££'^ m „(Q^ nrl e o, ©< 2 0 (-y. 


(23) 

(24) 


which is akin to Eshelby’s (1957) tensor for an ellipsoidal inclusion. 

Now 6Di„ rs = 0 if the y3th subvolume resides in the matrix, so that 

S£i% = 0 for M < p < N. (25) 

Because only M unknowns (associated with the subvolumes in the fiber) are involved in 
(23), we are left to solve 


M 


(26) 


A elf = &e° kl - l fStf rs AeJ 0 for a = 1, 2, A/. 

0=1 

When this relation is assembled columnwise for each fiber subvolume a, the solution 
can be obtained by Gaussian elimination. However, the square matrix which results from 
assembling these equations is of order 6M, and for a large number of subvolumes, M, may 
pose storage problems on the computer. Instead, we solve the equations by an iterative 
method. 


3.1. Magnification tensor 

If we single out the ath fiber subvolume on the right-hand side of (26), we can write 


M 


A elr = As* 0 , - rsx Ael a - £ /%f s Aejf, 


0 = 1 
0 ^ at 


or, on rearranging, 


a eir = ihm.n+rs^r'l asL - £ fs^ nn a £ ; 


T 3 


(27) 


(28) 


0 = i 

0 * a 


for a = 1 , 2, . . . , M, in which I k(mn denotes the fourth-rank identity tensor. This equation 
can now be solved by iteration in the form 

/ m 


SAejHx+j - [Ifcfmn + f a Skhnn\ (A £ mn Z f^mnrsi^, 


T3| \ 

rs i\ » 


( 29 ) 
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until the (A + l)th iterate differs insignificantly from the Ath ^terate. If _we take 

guess as (Ae™)' = Ae® , the first iterate (AcJ/l, yields the Rayle.gh-Born approximauon 

to the total strain increment. Continued iteration yields ’ higher-order Ryg ^ 

approximations which converge to left-hand side of 

dominant diagonal terms containing S mnrs in (27) and ta g 
the equation is required in order for the iterations to converge. 

We may write nm 

(AejHx+i = (30) 

where the operator (r*/„l 0 is given by 


u 


= ( hip, + /"W‘( W - Zf'S*, 


p = i 
(3 * a 


The operator (I^lx can be obtained recursively from the relation 


S<x ftta j - 1 


(Twrilx+l — Uklpq + f $klpq 


L fs% mn in„Jx 

0=1 , 

0*a ' 


(31) 


(32) 


which is obtained by combining (29) and (30). The magnification tensor MZ,„ for the ath 
fiber subvolume may then be written as 

M% lrs = l' m ITfc/nlx* ^ 

X 


and therefore the total strain increment in the ath subvolume is given by 


A elf = M\ 


klrs ^ b rs’ 


(34) 


Once the values of Aejf in the fiber (where 1 < a < M) are known, the values in he 
matrix (where M < a < N) can be found from (23). If further resolution is required, the 
value of AeJ r (r) at any point r in the unit cell can be found from (21). 

3.2. Rectangular subvolumes 

The iterative solution requires the evaluation of the tensor from (24). For 
isotropic constituents the tensor g klmn { Q 5lf mnrs may be written from (9) and (10) in the 
form 


SumJL Qt&m, 


X B - A' 


A m + 2 fi' 


SMi 


+ 


- 2 


R 0 - n m 


2/i ' 
.0 


+ S ksCrQ + KLtk + tlsCrtk) 


- v m \fr + /f” 


A m + 2//' 


•CrCsCfcCr. 


(35) 


where A 13 , n 0 and A m , /i m are the Lam6 constants for the /?th subvolume and the matrix, 
respectively, and 

& 


Ci = 


for / = 1,2 


(36) 


V(2t xn^Ttf + (27t« 2 /L 2 ) 2 

is a unit vector in the direction of the Fourier wave vector & = 2 nn-JU (no sum on i) 
The remaining factor required for assembling S ktrs is the Laue interference Integra 
product, Fo r the ath and /Jth rectangular subvolumes whose sides are ot 

length L°, L\ and L?, if in the x, , x 2 directions, we have 
Q° , (Z s )Q B (-fy = cos(4,(xf - xf) + Z 2 (X2 - xf)l 

sin(£.“^ 1 /2) si nfl"^^) sin(L^ t /2) sin(/. 2 ^2- / 2) (37) 

x (Ift,/2) (LIW 2) (Zfti/2) (iS« a /2) 

where xf, x 2 “ and xf, xf are the coordinates for the centers of the ath and /?th rectangular 
subvolumes. 
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4. NUMERICAL EXAMPLES 

Figure 2 shows the rransverse srress I.h 

periodic ceil, when a W/Cu occupies a volume fraction 

an overall stress of cr n - 1000 kPa. Z _ 395 GPa, v w = 0.28, 

/, 9/49 - 0 184 in the unh “"“ f u‘e 2 W presents a numerical tabulation of the 

£ cu - WGPz and c " ' h oftll e4 9 square subvolumes. These results are presented 
constant- valued stresses for each of the 4V square is a 

for those readers Altholgh Fig. 2(a) 

is°a correct graphical interpretation of the subv °[“™ ^ Jn the'fom ofcon” 
fields are uniform in — A^SSrtll* <han the 7 X 7 

<MS paper would lead to less interpolation error innicted by 

,h ' reXS'S^pproach-employing 

<* - -■ 9 — aK used “ 
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Fig. 2. Transverse stress concentration, <z„ for an appli 

Each uni, cell, with its 49 subvolumes, » J "j'^^S'Nunitad values for the 

The homogenized transverse Young s rnodulu in the center of the unit cell 

^ ssrs 

-- 5 -h rows 3, 4, 3. 
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calculate the stress variation throughout the tungsten fiber, whilst 40 subvolumes are used 
for the copper matrix. The stress concentration in the tungsten fiber varies from 1175 to 
1392 kPa, with a volume average of 1297 kPa. In the copper matrix the stress concentra- 
tion varies from a minimum of 682 kPa to a maximum of 1210 kPa. With 49 subvolumes 
in the unit cell, the overall homogenized transverse elastic modulus is calculated to be 
< e > = 156.3 GPa. It can be seen that the transverse stress in the square planform fiber 
forms a ridge/valley in the direction of the transverse stress. The average stresses in the 
ridges are 1340, 1392 and 1340 kPa, and the average stresses in the valley are 1179, 1 175 
and 1 179 kPa. This behavior can be noticed in a similar problem where a cuboidal inclu- 
sion in an infinite matrix suffers a uniform eigenstrain or transformation strain. The 
problem is outlined on p. 107 of Mura’s (1987) book “Micromechanics of Defects in 
Solids”. The ridge/valleys are present even when the fiber is isolated, and may be con- 
trasted to the case of an isolated fiber with circular (or ellipsoidal) cross-section which, by 
Eshelby’s (1957) analysis, would possess a uniform stress distribution within the fiber. 

A check on the transverse elastic modulus can be made by assuming that the unit cell 
is comprised of four subvolumes, with one subvolume in the fiber and three subvolumes 
in the matrix, as in Aboudi’s (1987) model. Each subvolume is assumed to be a spring so 
that the unit cell is assimilated to two springs connected in parallel. One spring in this 
parallel arrangement consists of a fiber and matrix spring in series, whilst the other spring 
in the parallel arrangement is a matrix spring. Taking account of the volume fraction of 
each spring allows the homogenized transverse modulus to be written in the form 

<£> = Ec(l - V + i . f (1 5 ( £c ,/ £w ) )) • <38) 

Assuming the volume fraction of the tungsten fiber to be/ = 9/49 gives < E ) = 149.3 GPa 
vs 156.3 GPa from the Fourier series calculation. Experiments carried out at NASA Lewis 
Research Center (Verrilli, 1988) have given values of the transverse elastic modulus of 
W/Cu composites as <£> = 136 ± 15 GPa at 9% volume fraction and <£> = 178 ± 15 
GPa at 40% volume fraction. An interpolation gives <£> = 148 ± 15 GPa at a volume 
fraction of/= 9/49 = 18.4%. 

Figure 3 shows the hydrostatic stress ( o u + o 22 + o 2i )/2> for the same transverse 
loading condition. High hydrostatic stresses occur at the fiber/matrix interface perpen- 
dicular to the direction of the loading axis. If the fiber had the same elastic moduli as the 
matrix and the unit cell was elastically homogeneous, the average hydrostatic stress in 
each subvolume would be 333 kPa. Tungsten/copper composites fail in thermomechani- 
cal fatigue tests through grain boundary cavitation in the copper matrix near the interface, 



1 2 3 4 5 6 7 


Fig. 3. Contour plot of the hydrostatic stress concentration, (cr M + a u + <7 33 )/3, for an applied 
transverse stress of a?, = 1000 kPa with a gradation of 50 kPa per contour. For a homogeneous 
material the hydrostatic stress would be 333 kPa. A fiber is located at the intersection of columns 
3, 4, 5 with rows 3, 4, 5. Each unit cell, with its 49 subvolumes, is embedded in a doubly-periodic 

array of identical cells. 
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Fig. 4. Contour plot of the transverse stress concentration, a Ut for an applied transverse stress of 
a,, = 1 000 kPa with a gradation of 200 kPa per contour. A fiber is located at the intersection of 
columns 3, 4, 5 with rows 3, 4, 5, except the (5, 3) node which is a void. Each unit cell, with its 
49 sub volumes, is embedded in a doubly-periodic array of identical cells. 



Fig. 5. Contour plot of the transverse stress concentration, <r n , for an applied transverse stress of 
Ou = 1000 kPa with a gradation of 300 kPa per contour. A fiber is located at the intersection of 
columns 3, 4, 5 with rows 3, 4, 5, except the (5, 3) and (5, 4) nodes which are voids. Each unit cell, 
with its 49 subvolumes, is embedded in a doubly-periodic array of identical cells. 



Fig 6 Contour plot of the transverse stress concentration, a n , for an applied transverse stress of 
o° = 1000 kPa with a gradation of 200 k Pa per contour. A fiber is located at the intersection of 
columns 3, 4 with rows 3, 4, 5, and a void is located at the intersection of column 5 with rows 3, 4 5. 
Each unit cell, with its 49 subvolumes, is embedded in a doubly-periodic array of identical cells. 
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and by tensile overload of the tungsten fibers (Kim et al . , 1989). High hydrostatic stresses 
near the interface may be important for cavitation growth in the copper matrix, which is 
an important creep rupture mechanism known to occur in copper. 

It is also of interest to examine the redistribution of stress when one or more of the 
nine sub volumes comprising the fiber is assumed to lose its load-carrying capacity, 
thereby becoming a void. One corner subvolume is assumed to be a void in Fig. 4. The 
transverse stress concentration in the rest of the fiber is now enhanced to compensate for 
the loss in the load-carrying capacity of the void subvolume; in particular, the peak 
stressed subvolume in Fig. 4 is 1800 kPa, compared with 1392 kPa in Fig. 2. Figures 5 and 
6 show the stress redistribution when two and then three fiber subvolumes at the 
Fiber/matrix interface lose their load-carrying capacity. In these figures, the peak stressed 
subvolumes are 2419 and 2040 kPa, respectively, with the latter occurring in the matrix, 
as compared with a peak stress of 1392 kPa in Fig. 2. When viewed in sequence, Figs 4, 
5 and 6 show how the transverse stress field could vary due to the growth of a fiber 
debond or a crack. A finer meshing would permit a more detailed study of such flaws. 

Although W/Cu composites have a strong thermodynamically-compatible bond at 
the interface (i.e. there is no interspecies diffusion), it was thought worthwhile to 
investigate the behavior under a transverse load when the interface is composed of a 
degraded material, or perhaps is coated with a compliant layer of a third material. 
Specifically, the central subvolume is assumed to be pure tungsten, and its eight nearest 
neighbors are assumed to have elastic moduli that are the average of those for tungsten 
and copper. The transverse stress concentration is shown in Fig. 7. The ridge/valley has 
disappeared, but this is perhaps because only one subvolume is considered for the tungsten. 
Also, the stress field of the unit cell is more uniform than that of Fig. 2, as expected. 



Fig. 7. Contour plot of the transverse stress concentration, ct, , , for an applied transverse stress of 
o° X{ = 1000 kPa with a gradation of 50kPa per contour. A fiber is located at the (4, 4) node. Its 
eight nearest neighbors (an interface or compliant layer) comprise a material whose elastic moduli 
are the average of those for tungsten and copper. Each unit cell, with its 49 subvolumes, is 
embedded in a doubly-periodtc array of identical cells. 


5. CONCLUSIONS 

A magnification tensor is derived using Fourier series techniques. This tensor trans- 
forms the far-field strain to the local strain of a constant-strained subvolume, which is 
considered to deform elastically, A set of these subvolumes can be configured by the 
analyst to construct a representation for the unit cell of a periodic composite. The Laue 
interference integral associated with the geometry of a rectangular subvolume is given. 
Numerical examples for a fibrous W/Cu composite are used to illustrate the theory. 
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